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INTRODUCTION 
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The  objective  o!  this  work  is  to  develop  fast,  accurate  techniques 
to  determine  the  characteristics  of  wire  antennas  over  plane  imperfect 
ground.  Since  Sommerfeld,  who  first  analyzed  the  problem  in  1908, 
numerous  papers  and  articles  have  been  written  by  several  different  authors. 
A good  bibliography  for  this  general  area  can  be  found  in  the  book  by 
Banos  [1].  Recently,  Miller  et  al  [2,3],  Lytle  and  Lager  [4,5]  developed 
user  oriented  computer  programs  for  analysis  of  the  electromagnetic  char- 
acteristics of  arbitrary  wire  configurations  over  imperfect  ground.  With 
the  aid  of  these  general  purpose  program:!  one  can  easily  solve  for  input 
impedances,  current  distributions,  and  field  patterns  of  wire  antennas 
operating  in  the  presence  of  a plane  imperfectly  conducting  earth.  This 
thesis  deals  with  some  useful  revisions  and  improvements  of  their  work. 


[1]  A.  Banos,  "Dipole  Radiation  in  the  Presence  of  a Conducting  Half- 
Space,"  Pergamon  Press,  New  York  1966. 

[2]  E.  K.  Miller,  A.  .1.  Pozzio,  G.  .1.  Burke  and  K.  S.  Selden,  "Analysis; 
of  Wire  Antennas  in  the  Presence  of  a Conducting  Half-Space:  Part  1. 
The  Vertical  Antenna  in  Free  Space,"  Canadian  .1.  of  Physics,  vol.  50, 
pp.  879-888,  1972. 

[3]  F..  K.  Miller,  A.  .1.  Pozzio,  G.  I.  Burke  and  F.  S.  Selden,  "Analysis 
of  Wire  Auitennas  in  the  Presence  of  Conducting  Half-Space:  Part  II. 

The  Horizontal  Antenna  in  Free  Space,"  Canadian  .1.  of  Physics, 
vol.  50,  pp.  2614-2627,  1972. 

[4]  D.  L.  Lager  and  R.  .1.  Lytle,  "Numerical  Evaluation  of  Sommerfeld 
Integrals,"  Lawrence  Livermore  Laboratory,"  Rept.  UCRL-51688, 

October  23,  1974. 

[5]  I).  L.  Lager  and  R.  I.  Lytle,  "Fortran  Subroutines  for  the  Numerical 
Evaluation  of  Sommerfeld  Integrals  Unt er  Aderemem,"  Lawrence  Livermore 
.aboratory,"  Rept.  UCRL-51821,  May  21,  1975. 
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The  major  problem  in  developing  a general  purpose  program  to  approach 
the  problem  at  hand  is  the  required  evaluation  of  certain  semi-infinite 
integrals  encountered  in  the  exact  Sonmerfeld  formulation.  In  this  work  an 
efficient  method  is  applied  to  integrate  the  semi-infinite  integrals  in  an 
essentially  exact  way.  The  method  applied  involves  the  use  of  orthogonal 
interpolatory  polynomials  of  highest  algebraic  accuracy.  Where  applicable, 
a modified  method  of  steepest  descent  has  been  applied  to  evaluate  these 
semi-infinite  Integrals  to  reduce  the  time  of  computation  without  significant 
loss  of  accuracy.  In  addition  a method  of  steepest  descent  has  been  used  to 
obtain  the  reflection  coefficient  method  as  the  leading  term  in  the  series 
resulting  from  the  integration  of  the  semi-infinite  integrals.  Since  the 
reflection  coefficient  method  is  very  easy  to  apply  it  is  useful  to  examine 
its  relative  accuracy. 

Finally  methods  are  presented  here  for  optimizing  certain  performance 
indices  of  arrays  of  arbitrarily  oriented  wire  antennas  operating  over 
imperfectly  conducting  ground.  The  indices  considered  include  directivity, 
maximum  power  gain,  quality  factor,  and  the  efficiency  Index  referred  to 
as  main  beam  radiation  efficiency.  Optimization  problems  both  with  and 
without  constraints  on  the  resulting  antenna  pattern  are  illustrated. 
Optimization  of  the  performance  indices  with  respect  to  leading  and  vari- 
ation in  the  antenna  structure  are  also  examined. 

The  contributions  of  the  work  of  this  thesis  include  the  following: 

1)  Application  of  interpolatory  quadrature  formulas  to  integrate  the 
semi-infinite  integrals  encountered  in  the  exact  Sommerfeld  formulation 
is  illustrated. 
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2)  Use  of  a modified  method  of  steepest  descent  to  evaluate  the  semi- 
infinite integrals  to  reduce  the  time  of  computation  without  signlf leant 
loss  of  accuracy  is  demonstrated. 

3)  The  commonly  used  reflection-coefficient  method  is  derived  in  detail 
and  Its  relative  accuracy  is  discussed. 

„ 4)  New,  user-oriented  computer  programs  for  treating  the  imperfect  ground 

problem  have  been  developed  and  described. 

5)  Formulas  and  results  convenient  for  application  to  optimization  and 
design  problems  have  been  presented. 


2.  VKRTICAL  DIPOLES  OVER  1 iE  HORIZONTAL  PLANE  SURFACE  OF  IMPERFECT  GROUND 

Analysis  of  radiatioi  from  parallel  elementary  vertical  dipoles  of 
the  electric  type,  situate*.  over  a horizontal  imperfect  ground  plane,  is 
discussed  in  this  chapter.  First,  the  exact  analysis  of  radiation  from 
the  dipoles  is  made  using  the  Sommerfeld  formulation.  The  semi-infinite 
integrals  encountered  in  this  formulation  are  evaluated  numerically 
in  an  accurate  way  using  orthogonal  polynomials.  For  field  points  far 
away  from  the  source  points  the  semi-infinite  integrals  are  evaluated 
using  a modified  saddle  point  method  to  reduce  the  computation  time  without 
any  significant  loss  of  accuracy.  The  reflection  coefficient  method  is  also 
derived  by  applying  a saddle  point  technique  to  the  semi-infinite  integrals. 

Finally,  a compariron  of  accuracy  is  made  between  the  reflection  coefficient 
method  and  the  exact  Sommerfeld  formulation. 

2.1.  Sommerfeld  Formulation 

An  elementary  dipole  of  moment  Idz  oriented  in  the  z-direction  is 
located  at  (x' , y',  z')  as  shown  in  the  Fig.  1.  The  dipole  is  situated 
over  an  imperfect  ground  plane  characterized  by  a complex  relative  dielec- 
tric constant  t It  is  possible  to  formulate  a solution  to  the  problem 
of  radiation  from  the  dipole  operating  in  the  presence  of  the  imperfect 
ground  in  terms  of  a single  Hertzian  vector  it  of  the  electric  type. 

A time  variation  of  exp(jiot)  is  assumed  throughout  the  analysis.  The 

Hertzian  vector  Q tt  in  this  case  satisfies  the  wave  equation  [6] 
z z 

1 

[6]  A.  Sommerfeld,  "Partial  Differential  Equations  in  Physics," 

Academic  Press,  New  York:  1964. 

.1 

I 
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(2.1) 


2 2 V 
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>eOel 
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(2.2) 
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“ Voe2 


(2.3) 


(2.4) 


The  electric  and  the  magnetic  field  vectors  are  derived  from  the  Hertzian 
vector  from 


and 


*►  -v  2-> 

E = V(V  • it)  + k it 


ft  = jue0er(V  x n) 


(2.5) 


(2.6) 


respectively.  At  the  interface  z = 0,  the  tangential  electric  and  magnetic 
field  components  must  be  continuous,  conditions  which  in  terms  of  the 
Hertzian  vector  components  can  be  written  as 


3 3ir„ 

lz  _ 2z 

E1  9y  E2  3y 


(2.7) 


9TTlz  ^2z 

E1  3x  E2  3x 
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Since  these  boundary  conditions  must  hold  at  z = 0 for  all  x and  y,  the 
x and  y dependence  of  the  fields  on  either  side  of  the  interface  must  be 
the  same.  Therefore, 


= ^2^22  [fr°m  (2.7)  and  (2.8)]  (2.11) 

3 7T 

= — r—  [from  (2.9)  and  (2.10)]  (2.12) 

d Z dz 


The  complete  solutions  for  the  Hertzian  vectors  satisfying  the  wave 
equations  (2.1)  and  (2.2)  and  the  boundary  conditions  (2.11)  and  (2.12) 
have  been  given  by  Sommerfeld  [6],  Wait  [7],  Miller  et.  al.  [2]  and 
others  [8-9].  The  solutions  are 


*lz  * P[exp(-jk1R1)/R1  + 


f jQ(xp)  c2Va2  - k2  - Ejl  A2-k2 


0 e2^2~kl  + el^  ^ " k 


exp{  -\A2-k2  (z+z'  ) } AdA]  (2.13) 


and 


'2,  = 2P 


„ I?  2 

f JQ(Ap)  exp  {*yA‘-k2  z 


->T  A"-kt  z'} 


.2  , 2 
C1 


AdA 


J 9 2 /~2 

0 A - k^  + t^A  - k 


(2.14) 


for 

Re[^  A2  - k2  „]  > 0 . Here 

IfZ 


[6]  A.  Sommerfeld,  "Partial  Differential  Equations  in  Physics,"  Academic 
Press,  New  York,  1964. 

[7]  R.  E.  Collin  and  F.  J.  Zucker,  "Antenna  Theory:  Part  2,  pp.  386-437, 
McGraw-Hill  Book  Co.,  New  York  1969. 

[8]  G.  Tyras,  "Radiation  and  Propagation  of  Electromagnetic  Waves," 

Academic  Press,  New  York,  1960. 

[9]  L.  M.  BrekhovBkikh,  "Waves  in  Layered  Media,"  Academic  Press,  New  York. 


ja)4TTe0e] 


(2.15) 


p -V  <x  - x')2  + (y  - y"*2 


(2.16) 


^ - ^J~p2  + (z  - z')2 


(2.17) 


The  primed  and  unprimed  co-ordinates  are  for  the  source  and  field  points, 

V 

respectively.  For  the  first  term  inside  the  brackets  can  be  inter- 

preted as  the  particular  solution  or  the  direct  contribution  from  the 
dipole  source  and  the  second  term  can  be  interpreted  as  the  complementary 

solution  or  a reflection  term  (reflection  from  the  imperfect  ground  plane). 

V 

Similarly  the  solution  for  can  interpreted  as  a partial  transmission 

V 

of  the  wave  from  medium  1 into  medium  2.  With  these  thoughts  in  mind,  ir_ 

lz 

can  be  split  up  into  two  terms 


V V direct  . V refl 

TT-  «=  7T_  + TT,  “ P * (g  + g ) 

lz  lz  lz  B0  B8 


(2.18) 


where 


^direct  = p . e*p(-Jk1R1)/R1  - P • gc 


(2.19) 


V refl  _ 
"lz  = P 


X2-k2  - X2-k2  JQ(Xp)  exp[-^A2-k2 (z+z' ) ] 


XdA  (2,20) 


0 F2V2-k2  + X2-k2 


r2  2 
\ - k^ 


P . g 


A physical  explanation  to  the  two  components  of  the  Hertz  potential  tt^ 
can  be  given.  The  first  one  can  be  explained  as  a spherical  wave  originating 


at  the  source  dipole.  This  term  is  easy  to  deal  with.  The  difficult  problem 

V refl  V refl 

lies  in  the  evaluation  of  n^z  . can  be  interpreted  as  a superposition 
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of  plane  waves  resulting  from  reflection  of  the  plane  waves  into  which  the 
original  spherical  wave  was  expanded.  This  arises  from  the  identity 


exp(-jk1R2)/R2 


Jq(ap)  expf-'V  A2  - k2(z  + z’)] 

l)  >2  - kj 


Re[ ) X2  - k2]  > 0 and 


R2  = V p2  + (z  + z’ )2  . 


The  term  under  the  integral  sign  of  (2.21)  can  be  recognized  as  plane 

wave  decomposition  of  the  original  spherical  wave  source.  Upon  reflection 

V refl 

of  the  plane  waves  from  the  dipole  source  as  expressed  in  ir.  , the 
amplitude  of  each  wave  must  be  multiplied  by  the  reflection  coefficient 
R(A).  The  complex  reflection  coefficient  R(A)  takes  into  account  the 
phase  change  as  the  wave  travels  from  the  source  (x',  y',  z')  to  the 
boundary  and  then  to  the  point  of  observation  (x,y,z).  The  reflection 
coefficient  R(A)  is  then  defined  as 


c2\j  A2  - k2  - A2  - k2 

£2\/x2-^  + C1\/>2-^ 


where  the  semi-infinite 


integral  over  R(A)  takes  into  account  all  the  possible  plane  waves.  As 
e2  ”,85  (2.20)  reduces  to  (2,21)  and  represents  a simple  spherical 

wave  originating  at  the  image  point.  This  physical  picture  will  later  be 
applied  in  the  derivation  of  the  reflection  coefficient  method  where  the 
effect  of  the  ground  plane  is  approximated  by  modifying  the  components  of 
the  plane  wave  decomposition  of  the  spherical  wave  originating  from  the 
image  of  the  source  dipole  but  multiplied  by  a specular  reflection  coef- 
ficient, R(8). 
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There  are  two  forms  oitt.  which  have  been  used  by  different  authors 

lz 

since  Sommerfeld  and  each  form  has  certain  distinct  advantages.  The  ad- 
vantages will  be  presented  in  a later  section,  but  the  two  forms  used  are 
presented  here. 


Vefl  - PlexpH^R^/R, 


lz 


- 2c 


X2-k2  JQ(Xp)  exp [-"J  X2-k2(z+z')] 


1 0 {x2-kl  c2^f?-k2  + 


XdX] 


Pl8l  “ 8sV] 


The  other  form  is 


V rpfl 

= P[-  exp(-jk1R2)/R2 


lz 


+ 2e, 


t JQ(Xp)  exp [ -)j  X2  - k2  (z  + z')] 


J r~2  T f~2 

0 e2  yX  -k^  + V ^ 


XdX 


(2.23) 


- P[-g,  + G ] (2.24) 

1 sv 

Both  forms  (2.23)  and  (2.24)  are  valid  and  their  equivalence  can  be 
observed  by  applying  (2.21). 

2.2.  Analysis  of  the  Reflected  Field 

The  field  in  the  upper  half  space  (z  > 0)  consists  of  the  direct 
field  from  the  source  dipole  situated  at  z m z' , its  image  situated  at 
z = - z' , and  the  correction  term  ggy  used  in  (2.23).  The  correction 
term  takes  into  account  the  nature  of  the  imperfect  ground,  because  as 


®sV  * °' 
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a)  Exact  Solution 


The  problem  of  determining  the  reflected  field  amounts  to  evaluating 
the  integral  ggy  since  gg  and  g^  are  easy  to  calculate.  To  do  this  a 
deformed  contour  similar  to  the  one  suggested  by  Miller  et.  al.  [2]  is 
chosen  and  is  indicated  in  Fig.  2.  The  deformation  of  the  contour  is 
possible  since  all  the  branch  points  and  the  poles  lie  on  the  second  and 
the  fourth  quadrants.  Hence  on  the  deformed  contour 


8sV  2^C1 


IQ(Xp)  exp[-j  yfx*  + k*(z  + z')]^\*'+k 


{> 


2 2 r~2  T 

^\J  X +k^  + X +k£ 


AdX 


+ 2e, 


f J^CBp)  exp - k2  (z  + z')] 


V 


2 2 
B -k2 


gdX 


2 2 
B -k1 


(2.25) 


where  B = >+  j5.  The  judicious  choice  of  5 as  the  limit  of  the  first 


integral  is  made,  because  then  the  poles  of  g at  + 

sv  — 


klk2 

■\£f 


have 


+ k 


negligible  effect  on  the  numerical  evaluation  of  the  integral.  The 
modified  Bessel  function  of  the  first  kind  in  the  first  integral  and  the 
Bessel  function  of  a complex  argument  in  the  second  integral  are  computed 
by  [10,  11].  Finally  the  two  integrals  are  computed  by  using  the  Gaussian 
interpolatory  quadrature  and  the  Gaussian  Laguerre  interpolatory  quadrature 
formulas,  respectively.  A description  of  the  method  is  discussed  in  Appendix  A. 


[JO]  T.  K.  Sarkar  and  J.  E.  Lewis,  "Accurate  Generation  of  Real  Order  and 
Argument  Bessel  and  Modified  Bessel  Functions,"  Proc.  IEE,  p.  34, 
January  1973. 

[11]  J.  E.  Lewis,  T.  K.  Sarkar  and  P.  D.  O' Kelley,  "Generation  of  Bessel 
Functions  of  Complex  Order  and  Argument,"  Electron.  Lett.,  vol.  7, 
No.  20,  pp.  615-616. 
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i mrinitmuaii 


It  was  because  of  these  methods  it  was  possible  to  evaluate  these  integrals 
in  reasonable  time. 

The  major  problem  of  this  technique  is  that  for  large  values  of  p,  the 
values  of  the  integrals  in  (2.25)  are  quite  large  and  they  have  to  be  subtracted 
from  each  other  to  yield  a relatively  small  value.  This  could  be  observed 
from  the  behavior  of  Ig(Ap),  since  it  is  a monotonically  increasing  func- 
tion of  p.  However  this  formulation  is  quite  useful  for  small  values  of  p. 


b)  Refl>xti?n  Coefficient  Method: 


For  large  values  of  p the  method  of  steepest  descent  is  applied  to 
analyze  the  fields  of  the  reflected  wave  at  distances  far  from  the  source. 
Under  this  circumstance,  it  is  convenient  to  use  an  infinite  integral 
representation  of  ggV  rather  than  the  semi-infinite  integral  representation. 
Use  of 


J (x)  - y[H(1)(x)  + H(2) (x) ] 
n L n n 

and 

HaW’>  - - H(2)<*) 

n n 


in  (2.23)  expresses  ggV  as 


(2.26) 


(2.27) 


8sV  Jr'l  j 


|-  V k2-A2  Hq2^(Ap)  exp[- j Jk2-  T2 (z  + z')] 
"/k2-A2  c?  /k2-A2  + z1  \f  k2-  A2 


AdA 


(2.28) 


where  C is  a contour  as  shown  by  Fig.  3.  ggV  is  now  a spectrum  of  plane 
waves  travelling  away  from  the  ground  plane  with  the  vertical  component 

I 2 2 

of  the  propagation  constant  as  <V  k^  - A . The  integral  in  (2.28)  contains 

n j 

double  valued  iunctions  „-A  . The  proper  sheet  of  the  double  valued 

functions  are  those  on  which  the  radiation  condition  [12]  is  satisfied, 
U2l  Ref.  76],  p.  189. 
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tail 


+ 0 as  z + ® for  a fixed  p. 


i*e.  g -*•  0 as  p + 00  for  a fixed  z,  and  g „ 

SV  sV 

Since 

H02)°P)  | Ap  | - S'  exp[-j(Ap-  ir/4)]  (2.29) 

the  convergence  of  (2.28)  is  assured  when  Im(A)  <_  0 as  p -+  ».  Convergence 
is  also  assured  for 


and  for 


If 


then 


and 


where 


T „ 2 ,2.1/2 

Im(k^  - A)  <0  as  z->°°. 


(ki  " x2)ll2TTt  ki  + Jki»  ki  > o,  kj  < o . 


(k2  - A2)172  = t'  + jr' 


k'  kV  - A' A" 
t"  = 


k'2  - kf - A’2  + A’’2 

T*  = f E (— ^ ) + (kjk'.J  - A'A"rr 


j^'2  i."2  ,,2  > it 2 

+ 1 1 j 1/2 


A = A'  + jA"  - Re(A)  + j Im(A) 


(2.30) 

(2.31) 

(2.32) 

(2.33) 

(2.34) 

(2.35) 


The  positive  sign  is  chosen  for  t'  since  t'  > 0 (from  (2.31)).  On  the 

path  of  integration  t"  < 0 and  if  k'^  ^ 0,  then  convergence  would  be 

assured  even  if  A"  = 0.  If  medium  1 is  lossless  k^  = 0 then  from  (2.33) 

A'  > 0 when  A"  > 0 and  A'  < 0 when  A"  < 0.  The  first  condition  A"  > 0 

conflicts  with  the  convergence  requirement  for  the  Hankel  function  '(Ap). 

This  problem  can  be  avoided  if  medium  1 is  assumed  to  be  lossy  (i.e.  , 

k"  f 0)  and  the  "lossless  case"  then  assumed  to  be  the  limiting  form  of  the 

/r'  \ 

expression  k'^  -*■  0.  Since  \Ap)  can  be  integrated  through  the  origin 
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( A * 0)  even  if  A"  ■ 0,  the  path  C of  Fig.  3 can  be  modified  to  the  path 
of  Fig.  4,  following  the  real  axis  from  -00  to  +°°. 

The  integral  in  (2.28)  can  be  simplified  by  making  the  following 
substitutions : 


sin  6 

(2.36) 

sin  0 

(2.37) 

* R2C0S  0 

(2.38) 

The  interpretation  of  the  angle  0 is  shown  in  Fig.  1.  Hence  by  the 
application  of  (2.29),  (2.36),  (2.37)  and  (2.38)  to  (2.28)  yields 


SsV*3 


exp[j{{-Tr/4  - cos  (B-0))]dg 


(2.39) 


where  k = anc*  ^ is  a path  in  the  complex  8 plane  as  shown  in 

Fig.  5. 

There  is  one  obvious  weakness  In  the  arguments  presented  to  derive 

(2.39),  namely,  that  there  are  points  on  the  path  for  which  the  argument 

(2) 

of  the  Hankel  function  (Xp)  used  in  (2.29),  is  not  large  and  may  even 
be  zero,  so  that  the  asymptotic  expansion  for  large  arguments  cannot  be 
used.  However  as  argued  by  Brekhovskikh  [9],  the  arguments  will  be  rigorous 
if  the  large  argument  approximation  is  used  only  after  the  path  of  inte- 
gration has  been  changed  to  the  path  of  steepest  descent  r^.  The  result  will 
then  be  the  same. 


Assuming  medium  1 to  be  lossless,  so  that  k^  is  real,  the  trans- 
formation A = k^  sin  B implies  (8  = 8*  + j8") 
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+ jx"  ■ k^(sin  g'  cosh  0"  + j cos  6*  sinh  g") 


X'  ■ k^sin  6'  cosh  g' 


(2.41) 


X"  = k^  cos  g*  sinh  8"  . 


(2.42) 


Hence  the  mapping  transforms  the  quadrants  of  the  X-plane  in  parallel 
strips  of  width  t/2  radians,  and  the  path  of  integration  from  X'  » -» 
to  * is  transformed  to  the  path  P^,  where  X"  = ImCk^sin  8)  = 0,  as  shown 

I 2 2 

in  Fig.  5.  The  requirement  Im(^  k^-X  ) < 0 amount  to  Im(k^  cos  8)  < 0 

on  the  path  of  integration,  or  for  k^  real  sin  8'  sinh  8"  > 0.  The  script 

1 U'  denotes  the  strips  of  the  8-plane  on  which  the  above  inequality  is 

satisfied  (upper  Riemann  sheet)  and  the  others  by  ' L'  (lower  Riemann  Sheet). 

The  path  then  totally  lies  on  U.  The  location  of  the  branch  points  at 

> = + k in  the  X-plane  are  now  transformed  into  sinB_  = f k in  the  8-plane 
1 1 


and  are  situated  at  + tt/2,  + 3tt /2  and  so  on.  The  branch  '*uts  Im(-^k^-Xi)  = 0 
are  now  transformed  to  Ira(k^cos  8)  = 0,  and  begin  at  the  branch  points 
P.g  . Since  Im(k^  cos  8)  = Im[k^sin(iT/2  + 6)],  these  branch  cuts  will  run 
parallel  to  the  path  1'^,  flm(k^,sin  8)  * X"  = 0]  but  shifted  by  + tt/2  along 
the  real  axis.  So  the  transformation  X = k^sing  has  transformed  the  upper 
and  lower  sheets  associated  with  the  branch  points  + k^  into  one  sheet 
where  certain  strips  on  the  sheet  belong  to  the  previous  upper  (U)  and 
lower  (L)  Riemann  sheet  on  the  X-plane. 

The  remaining  branch  points  X = + k are  transformed  to  sin  8 = + 

“ 2 ®2  ~ 

/ e2^ei  = i ^ which  has  solutions 


8 = j ln[+  j/tT  + A - e]. 


(2.43) 


The  branch  cuts  ImfCk2-*2)1/2]  = 0 are  transformed  into  Im[(e2  - e^sin  0)1/2]- 

1/2 

Im[(e  - sin  0)  ] = 0.  In  the  0-plane  there  are  then  two  Riemann  sheets 

1 /2 

connected  along  the  branch  cuts  Im[(e  - sin  B)  ] = 0 of  the  branch  points 
Bg  • Finally  the  poles  in  the  A-plane  are  now  given  by 

slnep'±ir?7  <2-44> 

2 

The  case  for  which  e=e'(l-j)  is  of  interest  since  this  is  a good  approxi- 
mation of  the  complex  refractive  index  of  both  earth  and  sea  water  over  a 
wide  range  of  frequencies.  The  displacement  current  is  negligible  in  com- 
parison to  the  conduction  current  if  the  frequencies  are  limited  below 

2 

about  1MHz.  For  a special  case  < = c'(l-j)  and  IUH  » 1,  the  locations 
of  the  branch  points  Bg  and  the  poles  Bp  can  be  approximated  by 


B 

B 

P 

P 


1,2 

2 

3,4 

2 


1,2 


2,3 


= + tt/4  + j ln(>/8e ' ) 

= + 3 tt/4  + j InC/ik5") 


[tt/2  + 1/(2  /P~)]  + 2 ^ 


- (tt/2  + — - — ) + — ^ — 
2/t1"  2/P~ 


and  they  are  illustrated  in  Fig.  5.  Out  of  the  possible  locations  of  the 

2 3 

branch  points  and  poles  only  , ?2  anc*  P3  are  situated  on  the  upper 

Riemann  sheet  of  branch  points  0 on  which  Im(k1cos  0)  < 0.  It  is  also 

B2  1 

important  to  note  that  none  of  the  poles  (P2>  P^)  are  situated  between  the 
original  path  of  integration  and  the  path  of  steepest  descent  Tq.  For 
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2 

this  particular  case  e = e'(l-j)  the  poles  Pj^  and  are  on  a line  having 

a slope  it/4  with  respect  to  the  6'  axis.  The  path  of  steepest  descent  also 

makes  an  angle  tt/4  and  it  increases  as  it  goes  up  in  plane.  Hence  the 

pole  will  never  be  intercepted  by  the  path  of  steepest  descent.  However 

when  the  path  of  steepest  descent  lies  in  close  proximity  of  the  pole 

P , special  precautions  must  be  taken  in  the  evaluation  of  the  integral 

of  (2.39).  The  pole  is  of  no  concern  since  it  lies  on  the  second 

Riemann  sheet  of  the  branch  point  B0  on  which  IrntJz  - sin^g)  > 0.  The 

B2 

presence  of  the  branch  point  should  ordinarily  be  taken  into  account 

while  deforming  the  path  to  the  path  of  steepest  descent  Tq.  The 
location  of  the  branch  point  is  such  that  its  imaginary  part  is  pro- 
portional to  In ( /8 e ' ) , which  for  large  e is  very  large.  Because  of 
this  reason,  as  explained  by  Tyras  [8],  the  contribution  along  the 
borders  of  the  branch  cut  would  be  a fast  decreasing  exponential  that 
can  be  neglected  in  comparison  to  the  contribution  from  the  saddle  point 
integration. 

For  0 < ,/2,  application  of  (B.2)  and  (B.4)  to  (2.39)  yields 


8sV 


„ 2 «-PHfci»2>  -J  c - sin26 
~ R 

4T  -2 


e cos  0 +\ 


sin  0 


_ 1 ^e(e-l)  [2c(e-l)-fEcos^0(3-cos^0)+cos  eVc-sin4Q(2c+sin^6)]  _ 1 ^ 

^R1R2  (c  - sin^0)^  [c  cos  0 + yj  e~-  sin2oT^  4sin^0 

(2.45) 


V ref  1 
*lz 


is  obtained  from  (2.23)  as, 
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I 


I 

I 


V ref  1 

TT, 

lz 


exp(-jk.R-)  . 

„ r J 1 2 ,£  cos  G 

Y . l 


-C 


2 

sin  0 


e cos  6 +>[c 

2 2 
r\  t r\  *-  r 


sin20 


JklR2 


{ 


e (e-1)  [ 2 e (e-1)  -4-  e cos^G  (3-cos^G)  + cos  6 (2e  + sin26fl/e  - sin20] 
2^3/2 


(c  - sin^G) J' “■  [c  cosG+VT^7]3 


V 


e - sin  6 


-}] 


2 sin20[e  cos0+VE  ~ sin20] 

For  very  low  dielectric  constant  this  method  would  he  inaccurate  as 

3 

tl  en  the  contribution  from  the  branch  cut  connecting  B2  to  (irt-j®)  would 
have  to  be  taken  into  account  and  the  second  terms  of  (2.45)  and  (2.46) 
would  be  different. 

The  first  term  of  (2.46)  which  has  a spherical  wave  representation 
of  the  original  source  but  now  emanating  from  its  image  can  then  be 
represented  as 


(2.46) 


M:fl~P  • rTM  exp(-jklR2)/R2 


(2.47) 


where  I'  can  be  recognized  as  the  specular  plane  wave  TM  reflection 
TM 

coefficient  [13],  given  by 


£ COS  0 


->T  E_- 


sin20 


TM 


£ COS  0 


+Ve- 


, 2„ 

sin  0 


(2.48) 


The  name  "Reflection  Coefficient  Method"  derives  from  (2.47)  since 
is  now  obtained  as  the  reflection  coefficient  times  the  potential  from 
the  image  of  the  source.  This  method  represents  a good  approximation, 
as  long  as  the  fields  are  computed  far  away  from  the  source  and  away 
from  the  imperfect  ground  plane  to  ensure  G < it/2.  The  far  fields  due 
to  the  reflected  ray  are  given  by  [14] 

[ 13]  E.  C.  Jordan  and  K.  G.  Balmain,  "Electromagnetic  Waves  and  Radiating 
Systems,"  pp.  628-654,  Prentice  Hall,  New  Jersey,  1968. 

[14]  R.  F.  Harrington,  "Time  Harmonic  Electromagnetic  Fields,"  McGraw- 
Hill,  New  York  1961,  p.  133. 
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Eref*  « jap.  sin  0 cos  0 P • 
p U lz 


and 


r.refl  4 . 2 V refl 

E = - 1o)|jn  sin  0 P ♦ tt. 
z J 0 lz 


or 


Erefl  = ju)pQ  sin  0 cos  0 P • exp^jk^) /R2 


= r_.  x (far  field  from  the  image) 
In 


(2.49) 


and 


E^ef  1 = - ju)P0  sin20  P • Tm  expf-jk^]/^ 


^TM*^ar  from  the  image) 


(2.50) 


The  direction  and  the  position  of  the  image  in  (2.49)  and  (2.50)  is 
given  by  considering  the  same  dipole  now  situated  over  a perfect 
ground  plane.  A physical  explanation  can  now  be  given  for  (2.49)  and 
(2.50).  A plane  of  incidence  is  now  defined  as  a plane  containing 
both  the  vertical  dipole  and  the  field  point  but  perpendicular  to  the 
ground  plane.  Then  the  contribution  of  the  reflected  ray  in  this  plane 

is  obtained  by  multiplying  E and  E (defined  in  this  plane)  by  the 

P 2 

specular  TM  reflection  coefficient.  This  is  different  from  the  exact 
solution  of  (2.20)  in  that  now  instead  of  infinitely  many  waves,  each 
multiplied  by  a reflection  coefficient  R(A),  all  the  plane  waves  are 
multiplied  by  the  specular  reflection  coefficient  R(0). 

When  the  conductivity  of  the  earth  is  large,  |e|  >>  1,  the  Hertz 
potential  in  medium  (1)  is  then  obtained  from  (2.18),  (2.19)  and  (2.45) 
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K Plexp(-jk1R1)/R1  + exp(-jk1R2)/R2  * 

[/c  cos  6-1  + 2e  ^ 1 j 3 + 

/T  COS  6 + 1 ^1^2  Jc  C08  0 +1 


}] 


Note  that  when  | c | 8oes  properly  Into  the  form  of  a field  due 
to  a vertical  electric  dipole  above  a perfectly  conducting  ground  plane. 
However  when  6 \ v/2 


V"lz  “ P[exp(-jk^R^)/Rj  - expC-jkjRj)/^  + 8e  2 expt-Jk^)*  ...)  (2.51) 

‘klR2 

The  sum  of  the  first  two  terms  may  be  smaller  than  the  third  term.  As  a 
matter  of  fact,  when  both  the  transmitter  and  the  receiver  are  on  the 
ground,  = R2  = p,  and  z = 0 = z’,  the  fields  will  be  solely  determined 
by  the  second  and  higher  order  terms.  The  reason  for  poor  convergence 
in  the  vicinity  of  8 * tt/2  is  that  the  effect  of  the  pole  at  tt/2  becomes 
important.  This  has  been  taken  into  account  in  the  formulation  of  the 
problem  in  the  next  section. 


c)  Fields  Near  the  Interface 


In  order  to  solve  for  the  fields  near  the  interface  a modified 
saddle  point  method  as  explained  in  Appendix  B is  applied  to  take  into 
account  the  pole  near  the  saddle  point.  The  expression  for  in 
(2.39)  has  a pole  Bp,  which  is  obtained  from  e cos  Bp  +^  e - sin20p  * 0. 
Hence 

sin  pp = VtVt  (2,52) 


and 

008  8p  ■ -frrr 


(2.53) 
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The  positive  sign  in  front  of  the  square  root  in  (2.53)  is  inadmissible 
because  then 

2 r 

Im(^e  - sin  8 ) = Im(- e cos  0 ) = - Im[  — * — J > 0 (since  e ■ e'  - je"). 

/e+1 

Hence  from  (B.5),  (B.6)  and  (B.10) 

2k  / 2~  I 2~ 

g5V=exp[-J»/4]  (^->  *-2— —f—==-  » 

2 e - 1 cos  6 - l.//e  + 1 


2 siri  (- 


0 + 8. 


exp [-  J^2.R2  C0S 
__ 

sin  ( — ) 


4nkjJ  i/2Vc  - sin20  Vc  - sln20-  c cos  0 
^2  - 1 cos  0 - l.//e  + 1 


where 


exp  I-  jk^  - W ] erfc(jW) 

+ co8  e + ^ sin  6 j 1/2 
/e  + 1 /"e  + 1 


2 ® 

- jk1R2  x 2 x sin  ( — ) 


jk  R [1  + COS-0  .jlliaJ., 
/e  + 1 /e  + 1 


(2.54) 


(2.55) 


V ref  1 

u,  can  now  be  obtained  after  g . is  obtained  from  (2.54).  This 
lz  sv 

method  then  can  be  applied  to  calculate  the  Hertz  potential  for  large 
p since  the  exact  solution  described  in  part  (a)  of  this  section  be- 
comes inaccurate.  This  method  is  used  to  computed  ^ for  a choice 
of  p > 0.5X.  The  primary  source  of  error  in  using  this  formula  for  a source 


(x'  = 0 , y ' = 0,  z'  = 0.25X)  and  a field  point  (x  = 0,  y = 0.5X,  z = 0.25X) 
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'■’Vue*-, , ...  ^ 


taffSJUT 


«MWW«m*w  w-’  »■  ■ 


(2) 

would  he  the  error  in  representing  Hq  (k^  p sin  6)  by  Its  large  argu- 
ment approximation  assuming  the  higher  order  terms  in  the  saddle  point 
method  are  negligible  compared  to  the  first  term.  For  this  case  the 
error  in  representing  by  its  large  argument  approximation 

is  one  percent  in  magnitude  and  -3“  in  phase.  But  ofter  we  would  be 
interested  in  antenna  arrays  closet'  to  the  ground  plane  than  0.25 A. 
Hence  the  error  would  always  be  less  than  the  upper  bound  presented. 

The  quantity  W presented  in  (2.55)  is  known  as  the  numerical 
distance  in  the  literature.  When  | r.  | is  large 


2 J^i^2  r—  9 

VTC=: [1  + (tt/2  - 0)  /T]2, 


(2.56) 


The  modulus  of  the  numerical  distance  W can  be  of  the  order  of  unity 

or  less,  in  spite  of  the  fact  that  k.^2  is  ^ar8e*  *s  c^us  evident 

that  the  condition  Iw2!  >>  1 can  turn  out  to  be  considerably  more 

rigid  than  the  condition  k^^  >>  1.  Different  approximations  used  for 

the  Hertz  potentials  mentioned  by  various  authors  [4,5]  in  literature 

are  generally  valid  for  |w|  >>  1 except  for  Banos  [1].  Somerfeld  [6] 

gives  a very  interesting  physical  interpretation  of  numerical  distance. 

It  is  the  difference  in  phase  shifts  k p and  k p of  a spatial  wave  and  a 

x s 

2 2 
k k 

2 1 2 

surface  wave  respectively,  where  k ■ — r — — . Sommerfeld  thus  defines  the 

8 k2+k2 

*1  *2 

numerical  distance  W1  to  be  = - j (^  - kg)p.  The  absolute  value  of  Wx  is 
small  compared  to  k^p.  In  fact  for  large  |e| 
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(2.57) 


k1p 


Is  I 


For  small  values  of  |W^|  the  spatial  wave  type  (due  to  and  g^)  pre- 
dominates In  the  expression  of  the  reception  intensity.  In  this  case, 
the  ground  constants  have  no  marked  influence.  For  larger  |w^| , 

the  relations  between  the  spatial  and  surface  waves  (due  to  g ) becomes 

S V 

important.  In  this  case  the  complex  ground  dielectric  constant  e has 
a marked  influence.  When  both  the  transmitter  and  the  receiver  are  on 
the  ground,  i.e.  0 css  it/ 2 then  | [ = | W.J  from  (2.55  and  2.57).  For 

n -rr/2  and  jw|  1,  the  disadvantage  associated  with  the  representa- 

•f  rpfl 

tion  of  ' by  (2.23)  is  that  the  second  and  higher  order  terms  in 

the  modified  saddle  point  evaluation  of  g^  cannot  be  neglected.  Even 
though  formulas  exist  for  generating  higher  order  terms,  the  procedure 
is  very  tedious  since  it  involves  taking  higher-order  derivatives  of  the 
integrand. 

The  reason  for  using  (2.24)  now  becomes  clear.  When  the  dipole 
and  the  observation  point  are  both  at  the  interface,  the  first  two 

V 

terms  gg  and  g^  of  in  (2.24)  tends  to  cancel  each  other  and  the 

field  is  solely  determined  by  the  integral  of  G , just  by  the  leading 

S V 

term  of  the  asymptotic  expansion  of  Ggy*  Substitution  of  (2.29),  (2.36), 

(2.37)  and  (2.38)  in  the  expression  G in  (2.24)  and  application  of  the 

sy 

modified  method  of  steepest  descent  as  of  (B.5)  and  (B.10)  yields. 
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sV 


g exp  (-j  tt/4) 


.2kl,ln  B,l/2 
irR2Bln  0 


cos  g 


c cos  8 + 


- sin% 


exp[-jk1R2  cos(g  - 0)  ] dB 

^>1T^l^l/2  cos  0 a/ e - sin20  - e cos  9 

R2  cos  0 - l./JTTT  e2  - 1 

exp[-jk  R - W2]  erfc(jW) 

x (2.58) 

{1  + c08  Q + /e  sin  6-|l/2 

/e  + 1 /c  + 1 


2 

where  W is  given  by  (2.55). 

Hence  for  |w|  > 1 or  k^R2  > °le|,  the  second  form  of  as  given 

by  (2.24)  is  used  to  compute  the  Hertz  Potential. 


2.3.  Evaluation  of  Mutual  Impedance  between  two  Parallel  Vertical 
Dipoles  Situated  over  a Plane  Imperfect  Ground 

In  this  section,  an  expression  is  obtained  for  the  mutual 
impedance  between  two  z-directed  current  elements  of  lengtht  l and  £ 

A D 

carrying  a current  distribution  1^  and  Ig,  situated  over  a horizontal 

imperfect  ground  plane.  The  corresponding  configuration  Is  presented 

in  Fig.  6.  The  electric  field  "£ . in  medium  1 from  the  current  element 

A 

IA  is  given  by  (2.5) 

E A = + V }Viriz  (replacing7*^  by  2.18) 

iwo  r 

' * ~ j 4 TA(8o  + 8s)dzi 
£A 

+ 1 i- 

ju)4ire  3z 


Vgo  + 8s>dzi 


(2.59) 
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7. 
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0.7 

h = X 

/fe^l 


Forbidden  region  for  the  reflection  coefficient 
method. 


medium  (1)  air 
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X 


medium  (2)  earth  y , c c,  k„ 

o o 2 


Fig.  6.  Two  parallel  z-directed  vertical  dipoles  of  lengths 

9 and  9 having  a current  distribution  of  I and  I„ 
An  Ad 

situated  over  a horizontal  earth. 


Since 


3gp 

o c 

3z  3z"! 


and 


SgB  8gs  . 

sT  'jiT  • and 

A 


assuming  the  current  I goes  to  zero  at  the  ends  of  the  open  wires, 

A 

3 3 

the  derivative  (r— ) operation  on  (g  + g ) can  be  transferred  to  r~ r 
3z  o s dz. 

A 

now  operating  on  I instead.  Hence 


t - 

A " Att 


2 IA(go  + 8s)  dZA 


j co4tt£ 


dIA 

4 <*„  - 8s>d2i 

A 


(2.60) 


fhe  mutual  impedance  Z„,  between  two  z-directed  current  elements  I.  and 

BA  A 

I is  expressed  as 
B 


BA 


^A  * ^B  dZB 


jwp( 


■b  dZB 


V8o  * gs)dZA 


j ojAttc 


XB  3Zg  { 


dz!  (go 
A 


gS)dZA,dZB 


(2.61) 


Transferring  the  derivative  operation  on  I in  the  second  term  and 

B 

assuming  I goes  to  zero  at  open  ends  of  the  wires 
B 

BA  4tt 


‘b  dZB  Vgo  + gs)dZA 


i r dIn  dIA 

j w4tte  J dzT  dzB  dzT  (go  “ gs^dZA 


0 o B 

B 


(2.62) 
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1 


If  a piecewise  sinusoidal  current  distribution  is  assumed  for  IA 
and  In  the  integrals  over  £.  and  £_can  be  performed  analytically.  But 
this  introduces  a third  order  pole  at  the  origin  of  the  semi-infinite 
integrals  defined  for  g . Hence  a four  pulse  approximation  to  the 
piecewise  linear  (triangle)  function  [15]  is  used.  In  this  case  also, 
the  double  integrals  in  (2.62)  can  be  evaluated  analytically,  leaving 
only  the  semi- inf inite  integral  to  be  evaluated  numerically. 

For  the  computation  of  the  self  impedance  Z , p in  (2.62)  is 
replaced  by  the  radius  of  the  wire  antenna. 

The  numerical  value  of  the  mutual  impedance  ZgA  or  the  self 

impedance  Z.A  lies  in  the  evaluation  of  g and  g . g is  exp(-jkR.)/R 
AA  O SO  11 

and  the  computation  is  straightforward.  The  different  techniques  discussed 
in  Section  2.2  are  utilized  to  compute  gg  either  in  the  form  (g^  - ggV) 
or  (-g^  + GgV) . In  the  reflection  coefficient  method  ggV  is  evaluated 
as  outlined  in  Section  2.2b.  In  the  "exact"  solution  ggV  is  integrated 
numerically  as  explained  in  2.2a  for  p < 0.5A,  and  evaluated  by  n modi- 
fied saddle  point  method  of  2.2c  for  p >_  0.5A.  For  p (2|e|/k^)A,  the 
(-g^  + G y)  representation  of  gg  is  used  in  the  computation.  User 
oriented  computer  programs  [16-19]  have  been  written  along  these  lines 


[15]  H.  H.  Chao  and  B.  J.  Strait,  "Computer  Programs  for  Radiation  and 
Scattering  by  Arbitrary  Configurations  of  Bent  Wires,"  Scientific 
Report  No.  7 on  Contract  No.  F19628-68-C-0180,  AFCRL-70-0374; 

Syracuse  University,  Syracuse,  New  York,  September  1970. 

[16]  B.  J.  Strait,  T.  K.  Sarkar  and  D.  C.  Kuo,  "Programs  for  Analysis  of 
Radiation  by  Linear  Arrays  of  Vertical  Wire  Antennas  Over  Imperfect 
Ground,"  Technical  Report  TR-74-1  on  Contract  No.  F19628-73-C-0047, 
AFCRL-TR-74-0042 ; Syracuse  University,  Syracuse,  NY:  January  1974. 

[17]  T.  K.  Sarkar  and  B.  J.  Strait,  "Analysis  of  Radiation  by  Wire  Antennas 
Over  Imperfect  Ground,"  Scientific  Rept.  No.  6 on  Contract  F19628-73-C- 
0047,  AFCRL-TR-75-0337 ; Syracuse  University,  Syracuse,  NY:  May  1975. 

[18]  T.  K.  Sarkar,  "Analysis  of  Radiation  by  Arrays  of  Vertical  Wire 
Antennas  Over  Imperfect  Ground  (Reflection-Coefficient  Method),"  IEEE 

Trans.  Ant,  and  Propagat.  vol.  AP-23,  September  1975,  p.  749. 

[19]  T.  K.  Sarkar,  "Analysis  of  Radiation  by  Arrays  of  Parallel  Vertical 
Wire  Antennas  Over  Plane  Imperfect  Ground  (Sommerfeld  Formulation) ," 
IEEE  Trans.  Ant,  and  Propagat.,  (to  be  published). 
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of  thought. 

Since  the  Reflection  Coefficient  Method  utilizes  a large  argu- 
ment approximation  of  the  Henkel  function,  maximum  error  in  this 
technique  occurs  in  computing  the  near  fields,  or  in  other  words  the  self 
impedance  Z . There  is  also  another  problem  associated  with  the  reflec- 

Aa 

tion  coefficient  method.  When  the  specular  angle  6 as  shown  in  Fig.  ]. 
approaches  90°,  the  saddle  point  0 approaches  the  pole  of  Fig.  5. 

Hence  the  reflection  coefficient  method  yields  a good  solution  when 
two  subsections  A and  B of  Fig.  6 are  away  from  the  ground  plane  and 
separated  by  an  appreciable  distance  yet  which  is  not  large.  It 
has  been  found  that  the  reflection  coefficient  method  yields  a result 
in  the  computation  of  impedance  elements  within  10%  of  the  "exact 
solution"  both  in  the  real  and  imaginary  parts  as  long  as  none  of  the 
elements  are  situated  below  a height  of  (0. 7//|  "e [ )1  from  the  ground 
plane  as  shown  in  Fig.  6.  This  also  ensures  that  when  two  elements  are 
as  much  as  1000A  apart  the  accuracy  of  both  real  and  imaginary  parts  of 
the  impedance  elements  computed  by  this  technique  are  still  within  10% 
of  the  "exact"  solution.  It  has  also  been  observed  that  the  restriction 
that  no  vertical  elements  can  lie  below  a height  of  (0.7//|e| )A  also 
ensures  an  accuracy  of  10%  in  the  impedance  elements  even  for  low  values 
of  dielectric  constant  (|cj^2)  of  the  ground  plane.  This  indicates 
that  so  long  as  this  restriction  on  the  heights  of  the  dipoles  is  main- 
tained the  contribution  from  the  branch  cut  of  is  not  too  important. 

2.4.  Conclusion 

Application  of  interpolator  quadrature  formulas  and  a modified 
method  of  steepest  descent  to  the  infinite  integrals  encountered  in 
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the  Sommerfeld  formulation  have  significantly  reduced  the  time  of 
computation  necessary  for  treating  the  imperfect  ground  problem 
without  significant  loss  of  accuracy  for  evaluating  these  integrals. 
When  a method  of  steepest  descent  is  applied  to  the  integrals  of  the 
exact  Sommerfeld  formulation,  the  reflection  coefficient  method  falls 
out  as  the  leading  term  in  the  series  under  certain  conditions.  The 
reflection  oefficient  method  yields  a result  accurate  to  within  10% 
of  the  exact  Sommerfeld  formulation  in  both  the  real  and  imaginary 
parts  of  the  impedance  elements  under  all  conditions  oi  the  ground 
so  long  as  they  are  (0.7//feJ)X  away  from  the  surface.  The  reflection 
coefficient  method  does  provide  a good  engineering  result  at  about  15% 
of  the  computation  time  taken  by  the  exact  Sommerfeld  formulation. 
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3.  HORIZONTAL  DIPOLES  OVER  A PLANE  IMPERFECT  GROUND 


In  this  section,  analysis  of  radiation  is  performed  for  par  Ilel 
horizontal  dipoles  situated  over  a plane  imperfect  ground.  As  before, 
for  near  fields  the  semi-infinite  integrals  obtained  in  the  Sommerfeld 
formulation  are  evaluated  numerically  using  orthogonal  interpolatory 
polynomials.  For  fields,  away  from  the  source  a modified  saddle-point 
method  is  applied  to  evaluate  the  same  integrals  to  reduce  the  time  of 
computation  without  any  significant  loss  of  accuracy.  The  reflection- 
coefficient  method  is  also  derived  as  a special  case  in  the  Sommerfeld 
formulation.  The  accuracy  of  the  reflection-coefficient  method  is 
compared  with  that  of  the  exact  Sommerfeld  formulation. 

3.1.  Sommerfeld  Formulation 

An  elementary  horizontal  electric  dipole  of  moment  Idx  and 
oriented  in  the  x-direction  is  located  at  a height  z'  above  the  hori- 
zontal imperfect  ground  plane,  as  shown  in  Fig.  7.  In  this  case  two 
components  of  the  Hertzian  vector,  n and  of  the  electric  type  are 
necessary  to  specify  the  fields  completely.  Assuming  a time  variation 
exp(jwt)  the  Hertz  potentials  satisfy  the  following  wave  equations  [6] 

(V2  + k2)^  = - R S(x-x')  6(y-y')<5(z-z')  (3.1) 

and 

(V2  + k2)r2  = 0 (3.2) 

The  electric  and  the  magnetic  field  vectors  are  given  by  (2.5)  and 
(2.6)  respectively. 
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A horizontal  dipole  over  an  imperfect  ground  plane. 


From  the  continuity  of  E^  at  the  interface  z * 0 it  follows  that 


3 rr_  3it, 

lx  lz 

3x  3z 


3 TT 


2x 


3 TT 


2z 


3x 


3z 


(3.3) 


Similarly  from  the  continuity  of  and  (3.3)  it  follows  that 


kl  ’lx  * k2  "2x 


(3.4) 


From  the  continuity  of  the  magnetic  field  component  at  the  inter- 
face z = 0, 


3tt 


lz 


3 TT 


2z 
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(3.5) 


and  from  the  continuity  of  it  follows  that 


,3’lx 
E1  ( 3z 


3 TT 
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3x 


,^2x 

e2  ( 3z 


3 TT 
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3x 


(3.6) 


Since  (3.5)  holds  for  all  values  of  y, 


C1  "lz  = C2  V2z 


and  because  of  (3.7),  (3.6)  reduces  to 


Sir,  3tt0 

lx  2x 

“ e. 


'1  3z 


'2  3z 


(3.7) 


(3.8) 


After  the  components  ir  have  been  solved  for  by  using  the  wave  equa- 
tions (3.1)  and  (3.2)  and  the  boundary  conditions  (3.4)  and  (3.8),  the 
components  can  be  obtained  from  (3.3)  and  (3.7).  The  complete  solu- 
tion for  the  Hertzian  vectors  has  been  given  by  Sommerfeld  [6],  Wait  [20] 
Miller  et.  al.  [3]  and  others,  and  can  be  written  as  follows: 

[20]  T.  R.  Wait,  "The  Electromagnetic  Fields  of  a Horizontal  Dipole  in 
the  Presence  of  a Conducting  Half-Space,"  Canadian  J.  of  Phys., 
vol , 39,  pp.  1017-1028,  1961. 
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v lx  " p' (exp(-jk1R1)/R1 
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+ | exp[-yX  -k^z+z')] 
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AdX] 
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JQ(Xp)expt'|/ X2-k2  z-V  X2-k2  z']XdX 


(3.12) 


for  Re[ >/x2  - k2  2]  > 0.  Here  p and  R^  are  given  by  (2.16)  and  (2.17), 


and 


P' 


IAx 


ja^ire^ 


(3.13) 


V H 

Like  tt^z  in  (2.18),  can  similarly  be  split  into  two  parts 


H H direct  . H refl 

IT,  = TT,  + 7T, 
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where 


^direct  = p,  . exp(-jk1R1)/R1  = P'  • g£ 
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(3.15) 


The  superscript  H represents  the  Hertz  potentials  due  to  a horizontal 
dipole,  whereas  V represents  the  potentials  due  to  a vertical  dipole. 
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- - - - 
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and  R2  is  given  by  (2.22).  Here  H7T^rect  represents  the  particular 
solution  of  the  wave  equation  (3.1)  or  the  direct  spherical  wave  con- 
tribution from  the  dipole  source  and  can  ^nterPretec^  as  t'ie 

complementary  solution  or  a reflection  term  (reflection  from  the  medium). 
It  represents  a superposition  of  plane  waves,  resulting  from  the  reflec- 
tion of  the  plane  waves  into  which  the  original  spherical  wave  from 
the  horizontal  dipole  source  was  expanded.  Upon  reflection,  the 
amplitude  of  each  plane  wave  must  be  multiplied  by  a reflection 
coefficient  V(A)  quite  different  from  R(A)  as  in  the  pre- 
vious section  because  of  a different  orientation.  The  complex  reflec- 
tion coefficient  V(X)  takes  into  account  the  phase  change  as  the  wave 
travels  from  the  source  (x‘ ,y',z')  to  the  boundary  and  then  to  the 
point  of  observation  (x,y,z) . The  reflection  coefficient  V(A)  is  then 


defined  as 


+ V*2-k 


As  before, 


2x 


can  be  interpreted  as  a 


partial  transmission  of  the  wave  from  medium  1 into  medium  2.  Unlike 

the  elementary  vertical  dipole,  the  horizontal  dipole  over  an  imperfect 

ground  plane  needs  two  Hertz  potentials  to  completely  specify  the  fields. 

|| 

The  primary  excitation  tt^  in  this  case  gives  rise  to  a secondary 

excitation  ir  whose  contribution  is  highly  dependent  on  the  complex  dielec- 

11 

trie  constant  of  the  ground  (e)  This  is  because  as  e °°,  ir,  0. 

lz 

|| 

Also,  as  e -*■  1,  0.  To  get  a physical  interpretation  of  the  nature 

H 

of  the  fields  due  to  is  useful  to  consider  the  following  form 
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(3.16) 


X2dX 


where  $ = (x,p) 
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A principal  distinction  of  the  horizontal  dipole  as  compared  to  the 
vertical  dipole  is  its  characteristic  radiation  implied  by  the  factor 
cos  <f>.  So 


2 H 

*1*  niz 
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J rr,  2 , n/TF“ 

0 e -y  X -k^  4-  yX  -k 


XdX  (3.17) 


where  e = e2^ei 


The  first  term  in  the  above  expression  can  be  interpreted  as  fields 
originating  from  a vertical  octupole  of  strength  (1  - ^-)  situated  at 
the  image  of  the  original  horizontal  dipole.  The  second  term  consists 
of  infinite  plane  waves  reflected  from  the  ground  plane  and  each  multi- 
plied by  a different  reflection  coefficient  R(X).  It  represents  a 
superposition  of  plane  waves  resulting  from  the  reflection  of  plane 
waves  into  which  a spherical  wave  from  a vertical  octupolar  source  of 
strength  (1  + was  expanded.  This  representation  is  in  no  way  unique. 
Interestingly,  for  large  values  of  t and  when  the  horizontal  dipole 
source  is  close  to  the  ground  the  octupolar  source  which  gives  rise  to 
reduces  to  a quadrupolar  source  as  illustrated. 
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(3.18) 
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(3.20) 


Thus,  for  the  source  close  to  the  ground  plane  and  for  large  e,  the 

horizontal  dipole  appears  as  a quadrupole.  Also,  tt_  can  be  in- 

Zz 

terpreted  as  representing  the  waves  induced  in  medium  2 by  the 

secondary  excitations  generated  in  medium  1 and  represented  as 
H 

IT  . 

lz 


3.2.  Evaluation  of  Mutual  Impedance  Between  Two  Parallel 
Horizontal  Dipoles 

In  this  section  an  expression  is  obtained  for  the  mutual 
impedance  between  two  x-directed  t urrent  elements  of  lengths  1 and  Jl 

A D 

carrying  current  distributions  1^  and  Ig,  respectively,  and  situated  over 
a horizontal  ground  plane.  The  corresponding  configuration  is  presented  in 
Fig.  8.  The  electric  field  in  medium  1 due  to  the  current  element  1^ 
is  given  by 


EA  * Ik!  + 7 7 


3tt 


2 . . 

k1  (».  + TT  ) + ( 

1 lx  lz  9x 
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8ir 
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8z 


(3.21) 


lx  lz 

For  an  elemental  dipole  of  moment  IAx,  the  expresf^on  (— — + — — )can 

dX  a Z 

be  written  in  a simpler  form.  By  using  (3.11),  (3.14)  and  (3.15), 


V 


Fig.  8.  TVo  parallel  x-directed  horizontal  dipoles  of  lengths 
2.  and  2 having  a current  distribution  of  I and  In 
situated  over  a plane  earth. 
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3x  + 3z  " P 3x  [g0  " 81  + (gsH  + ?z 


p*  . — [o  _ g + 2k2  — 

3x  l80  81  1 


JQ(Xp)  exp[-y  X2-k^  (z  + z')] 


0 k l^2-^  + k.2  y X2-k2 


= P'  * 3^  Ig0  ‘ 81  + GsV/c]  {fr°m  2>24) 


^ P 3x  tg0  “ gl  + g3 ^ 


Hence  the  x-directed  field  for  a current  distribution  I is  given  by 

A 


(3.21) 


x J**0 

EA VT  TA(g0  ' 81  + gsH)dXi 


+ -i  I_i_ 

3x  j47rwe 


rA^  f(80  ' 81  + 83)idxi 


since 


3gQ  3gx  3gx  3g3  3g3 

3x'  ’ 3x  3x'  anc^  3x  3x' 

A A A 


and  assuming  the  current  I goes  to  zero  at  the  ends  of  the  open  wires, 

A 

3 3 

the  derivative  — can  be  transferred  to  - — r operating  on  I.  so  that 

dX  dX.  A 

A 


jwu 

pX  _ o 

A 4tt 


V80  - *1  + 8sH)dxi 


+ h !i4=r  J s?  <go  ■ 8x + 83)  dxA^ 

° A 


The  mutual  impedance  Z„.  between  the  current  element  I.  and  another 

BA  A 


current  element  I oriented  in  the  same  direction  (ft)  is 

D 
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ea  'b 


jajpc 

4tt 


| IBdXB 


V8!)  ‘ 81  + 8sH^dXA 


nine  Id  dx°  ^ 

ft 


j4™e  J B B 3x 

0 l i 

B A 


dlf  (80  - 8l  + 83M*i 

A 


(3.25) 


By  transferring  the  derivative  to  I and  assuming  I goes  to  zero  at 

B B 

the  end  of  the  wires, 


'BA  4ft 


XB  dXB 


rA*80 


81  + 6sH)dxA 


ai. 


j4fra)e 


% 


dIA 

4feo 

A 


81  + «3)dxI 


(3.26) 


computation  of  mutual  impedance  between  two  parallel  horizontal  dipoles 
over  an  imperfect  ground  plane  then  amounts  to  evaluating  the  integrals 
of  g and  g_.  This  is  dealt  with  in  the  next  section. 

SH  j 


3.3.  Numerical  Evaluation  of  the  Integrals 

In  this  section  the  integrals  of  ggH  and  are  evaluated  in  an 
essentially  exact  manner.  Where  applicable  the  method  of  steepest 
descent  has  been  applied  to  reduce  the  computation  time  without  any 
significant  reduction  or  loss  in  accuracy. 

a)  Exact  Solution.  The  problem  of  determining  the  mutual  impedance 
amounts  to  evaluating  the  semi-infinite  integrals  gg^  and  g^in 
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an  essentially  exact  manner.  Since  the  Integrals  are  very  similar  to 
those  of  (2.23)  and  (2.24),  the  Integration  scheme  as  explained  in 
section  2.2a,  [3,  17)  is  used. 

One  of  the  major  problems  associates  with  this  technique  is  that 
for  large  values  of  p there  are  large  numerical  errors,  as  explained  in 
2.2a.  The  other  problem  is  that  when  the  current  elements  are  below  a 
height  of  0.03X  from  the  ground  plane  the  integrals  ggH  and  g^  in  (3.11) 
and  (3.15)  become  oscillatory  and  hence  the  quadrature  formulas  explained 
in  Appendix  A fail.  Therefore  it  is  not  possible  to  analyze  antennas 
situated  at  distances  less  than  0.03X  from  the  ground  plane  using  this 
procedure. 


b)  Modified  Method  of  Steepest  Descent:  For  large  p a modified  method 

of  steepest  descent  can  be  applied  to  the  Integral  g^  (because  of  a pole 
in  the  integrand)  and  a method  of  steepest  descent  can  be  applied  to  the 
Integral  ggH  to  reduce  the  computation  time  without  significant  loss  in 
accuracy.  Making  the  following  substitutions,  as  suggested  by  (2.26), 
(2.27),  (2.36),  (2.37)  and  (2.38),  the  integrals  can  be  written  as 


5sH 


2k  sin  8 1/2 

( — ± \ 

vttR  sin  6' 


cos  8 


cos  8 


exp[ j{-n/4  - 


sin^8 


kjR2  cos(8  - 0) })  d8 


(3.27) 


.2kl*ln  6 1/2 

^irR^in  0^ 


cos  8 


e cos  8 + 


2 

sin  8 


exp[  j(-ir/4  - 


kjR2  cos(B  - 0) )]  d0 


(3.28) 


where  1^  is  a path  in  the  complex  8-plane  as  shown  in  Fig.  5.  The  large  argu- 
ment approximations  for  the  Hankel  functions  have  been  justified  [9].  The  above 
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integrals  are  then  evaluated  by  using  (B.4)  for  g^  (since  the  expres- 
sion of  g^H  does  not  have  any  pole  near  the  saddle  point  6 ■ 6)  for 
gj,  respectively.  Hence, 


2 cos  6 


expC-jk^) 


ssH 


cos  0 +]j  c - sin2 0 


(3.29) 


[taking  only  the  spherical-wave  part  or  the  first  term  of  the  saddle- 
point  method]  and 


cos  6 


- 

g3  ( R,  cos  0- 

2 A/c+1 

where  W2  is  given  by  (2.55) 


sin  0 - e cos  0 


e2  - 1 


exp[-jk1R2~W2]erfc(jW) 

+ cos  0 { /c  sin  6 ^1/2 
/e+1  /e+1 

(3.30) 


An  user-oriented  computer  program  has  been  written  along  these 
lines  of  analysis  and  is  available  [17,21]. 


3.4.  Reflection  - Coefficient  Method 

For  field  points  far  away  from  the  source  and  away  from  the 
ground  plane  the  method  of  steepest  descent  is  applied  to  evaluate  the 
integrals  gg^  and  gy  For  both  source  and  field  points  near  the  ground 
plane,  the  effect  of  a pole  which  may  exist  near  the  vicinity  of  a 
saddle  point  (in  g^)  is  neglected.  Making  the  following  transformations 
as  suggested  by  (2.26),  (2.27),  (2.36),  (2.37)  and  (2.38),  and 

can  be  expressed  as 


[21]  T.  K.  Sarkar,  "Analysis  of  Radiation  by  Arrays  of  Parallel  Hori- 
zontal Wire  Antennas  Over  Plane  Imperfect  Ground  (Sommerfeld 
Formulation),"  IEEE  Trans.  Ant,  and  Propagat. , (to  be  published). 
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H refl  P' 

IT  * — - 

lx  2 


r 2k  sin  B 

(~'R  8^~n~  q-)  ' exp[j(-Ti/A  - k^RjCosCB  - 0)) 


cos  g — ~\/g  - sin^B 
cos  B +'\[c-  sinB 


dB 


IT-  ■ P 'cos  $ 
lz 


exp  { J {-  ./4-kjRjCoe  ( 6-6) } ) (,R_sln0 


2k,sinB  , /-  „ ,/"  , 2- 

1 v 1/2  cos  B -yt-sin  B 


)■ 


e ^08  B +>/e- sin2B 


X cos  BdB 


since 


h[2)(xp) 


Xp 


^TTXp’^1^2  exPM(*P  “ 3ff/4)] 


The  contour  of  integration  is  defined  in  Fig.  5.  Application  of 
the  saddle-point  method  as  given  by  (B.4)  to  (3.31)  and  (3.32)  yields 


H refl  cos 

TT,  “ P 

lx 


9 -Ve-Sin2e  “P<-JV2> 


"2 


cos  8 + ve-sin  9 


H xi*  * cos  0 -Ve.  - sin20  exp^ 

it,  * 2P  cos  d>  sin  0 cos  0 

lz  T , — R- 


/ 2~  "2 
e cos  0 +V  e - sin  0 


The  expressions  for  the  far  field  due  to  the  reflected  rays  only,  are 
given  by  [14], 

„refl  , . rH  refl  „ , H . „ . 

= - imp  cos  0 l rr,  cos  0 cos  6 - w,  sin  0J 
P J o lx  T lz 


..refl  . . - rH  refl  n . H . 

E = juiy  sin  0 [ IT-  cos  0 cos  ij>  - tt  sin  0] 
z O XX  xz 


..refl  H refl 

E,  - jwy  sin  it, 

4)  J o lx 


(3.31) 

sin  B 

(3.32) 

(3.33) 

(3.34) 

(3.35) 

(3.36) 

(3.37) 


(3.38) 


thi  i defined  plane)  by  the  specular  plane  wave  TM  reflection  coef- 
ficient and  (-  field  perpendicular  to  this  defined  plane)  by  the 
negative  of  the  specular  plane-wave  TE  reflection  coefficient.  The 

fields  E , E;  and  E used  for  this  case  are  assumed  not  to  originate 
p z 6 

from  the  source  dipole,  but  from  the  image  of  the  source  dipole. 

The  mutual  impedance  between  two  parallel  horizontal  electric 
dipoles  over  an  imperfect  ground  plane  can  be  computed  as  indicated. 
An  user-oriented  computer  program  is  available  using  the  reflection- 
coefficient  method  for  treating  problems  of  this  type  [22,23]. 


3.5.  Comparison  of  Accuracy  Between  the  Reflection  Coefficient 
Method  and  the  Exact  Sommerfeld  Formulation 

In  the  case  of  the  horizontal  dipoles  only  tt^z  has  a singularity 
and  not  tt^.  But  in  (3.16)  has  a cos  <J>  variation  which  indicates 

that  the  effect  of  the  singularity  would  be  maximum  when  the  field  due 
to  the  dipole  is  computed  in  a plane  containing  the  dipole  and  perpen- 
dicular to  the  ground  plane.  Wnen  the  field  is  computed  along  a plane 
perpendicular  to  the  axis  of  the  dipole,  the  effect  of  the  singularity 
is  zero.  Also,  from  (3.22),  the  contribution  from  the  singularity  is 
less  for  the  corresponding  vertical  dipole  by  a factor  of  1/e.  Hence 


[22]  T.  K.  Sarkar  and  B.  J.  Strait,  "Analysis  of  Radiation  by  Arrays  of 
Parallel  Horizontal  Wire  Antennas  Over  Imperfect  Ground,"  Scientific 
Report  No.  5 on  Contract  F19628-73-C-0047,  AFCRL-TR-74-0538,  Syracuse 
University,  Syracuse,  New  York:  September  1974. 

[23]  T.  K.  Sarkar,  "Analyi.  * of  Radiation  by  Arrays  of  Horizontal  Wire 
Antennas  Over  Imperfec..  Ground,"  (Reflection  Coefficient  Method)," 
IEEE  Trans.  Ant,  and  Propagat.,  (to  be  published). 
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Let  the  specular  TE  and  TM  reflection  coefficients  be  denoted  by 

T and  r , respectively,  and  defined  by 
It  TM 


cos  0 - V e - sin  6 


cos  0 + V £ “ sin  0 


(3.39) 


and  by  (2. 48). 
Since, 


H refl 

it,  cos  0 cos 
lx 


- tt.  sin  0 
lz 


= - P*  cos  0 cos  <j)  r^expC-jk^)/^ 


the  far  field  expressions  in  (3.36)  - (3.38)  reduce  to 


E^ef  1 = jwuo  cos20  cos2<f>rTMP'  exp(-jk1R2)/R2 


(3.40) 


(3.41) 


I’tm  • (far  Ep  field  from  the  image  due  to  a 
perfectly  conducting  ground  plane) 


- jwyo  sin  0 cos  0 cos  <J>  P'  exp(-jk1R2) /R2 


r_„  • (far  E field  from  the  image  due  to  a 

TM  Z 

perfectly  conducting  ground  plane) 


rpf  1 

E^  = ja)Po  sin  tp  rTE  P'  exp(-jk1R2)/R2 


= - T™  • (far  E,  field  from  the  image  due  to  a 

lcj  <p 

perfectly  conducting  ground  plane).  (3.42) 

A physical  interpretation  can  be  given  for  these  formulas.  A 
plane  of  incidence  can  be  defined  as  a plane  pasning  through  the  mid- 
point of  the  horizontal  dipole  and  the  field  point  in  question  and  that 
is  also  perpendicular  to  the  ground  plane.  Then  the  contributions  from 

the  reflected  ray  can  be  obtained  by  multiplying  E and  E (-  fields  in 

P z 
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it  is  expected  that  the  reflection-coefficient  method  would  give  a better 
approximation  to  the  exact  Somme-feld  formulation  than  for  the  vertical 
dipole  case.  Also,  as  before,  for  the  reflection-coefficient  method  the 
maximum  error  occurs  in  the  computation  of  the  near  fields  and  hence,  the 
self  impedances.  It  has  been  found  that  for  the  worst  case  of  two  coplanar 
horizontal  dipoles  (for  $ * 0)  the  reflection-coefficient  method  yields  a 
result  within  10%  of  the  exact  Sommerfeld  formulation  both  in  the  real  and 
imaginary  parts  of  impedance  elements,  as  long  as  the  dipoles  are  away 
by  at  least  (0. 25//[e"j")  X from  the  ground  plane.  For  the  reflection  coef- 
ficient method,  there  is  still  another  restriction  on  the  distances  between 
two  horizontal  dipoles.  This  is  because  as  the  distance  gets  larger  the 
saddle  point  0 approaches  the  pole  at  tt/2  [in  g^  (3.22)].  It  has  been 
observed  that  the  horizontal  dipoles  should  lie  within  a cone  whose  semi- 
vertical angle  is  70°  in  order  that  the  reflection-coefficient  method 
yield  results  within  10%  of  the  exact  Sommerfeld  formulation  both  in  the 
real  and  imaginary  parts  of  the  impedance  elements.  The  apex  of  the  cone 
is  the  image  of  the  middle  point  of  the  source  dipole.  This  is  illustrated 
in  Fig.  9.  This  implies  that  the  maximum  lengths  of  the  horizontal  dipoles 
for  which  10%  accuracy  can  be  obtained  in  the  evaluation  of  the  impedance 
elements  by  the  reflection-coefficient  method  is 

L = 2h  tan  70°;=5.5h  (3.43) 

critical 

as  indicated  in  Fig.  9. 

Also  for  large  e and  large  distances  away  from  the  source 

and  Hir,  in  (3.18)  and  (3.19)  can  be  rewritten  as 
lz 
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IT,  ■ P 
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k^/F 


z 3 exp(-jkiR) 
R 9R  R 


(3.44) 


H 


, -2j  p 3 

lz  k ^ R 3R  R 


(3.45) 


H H 

Hence  from  (3.44)  and  (3.45)~  and  tt^z  have  the  ratio  z/p  which  is 
very  small  at  great  distances  from  the  transmitter.  The  field  of  trans- 
mission of  a horizontal  antenna  has  the  same  character  as  the  field  of 
transmission  of  a vertical  antenna  except  for  the  $ dependence.  It  is 
then  expected  that  the  restriction  on  the  length  of  a horizontal  an- 
tenna can  be  removed  if  it  is  situated  at  least  (0.7//["Ff)A  away  from 
the  ground  plane  (as  from  2.3).  For  the  worst  case,  <f>  *»  0,  it  has  been 
found  that  the  reflection-coefficient  method  yields  a solution  within 
10%  of  the  exact  Sommerfeld  formulation  both  in  the  real  and  imaginary 
parts  of  impedance  elements  even  when  two  horizontal  currents  are  as 
much  as  1000A  apart. 


3.6.  Conclusion 

Application  of  interpolatory  quadrature  formulas  and  a modified 
method  of  steepest  descent  to  the  infinite  integrals  encountered  in  the 
Somerfeld  formulation  have  significantly  reduced  the  time  of  computation 
without  appreciable  loss  of  accuracy  for  evaluating  these  integrals.  When 
a method  of  steepest  descent  is  applied  to  the  integrals  of  the  exact 
Sommerfeld  formulation  the  reflection-coefficient  method  falls  out  as  the  lead 
ing  term  in  the  series  under  certain  conditions.  The  reflection-coefficient 
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method  yields  a result  within  102  of  the  exact  Sonmerfeld  formulation 
in  both  the  real  and  imaginary  parts  of  the  impedance  elements  pro- 
vided the  antennas  are  no  longer  than  5.5h  where  h,  the  height  of  the 
antennas  above  the  ground  plane,  can  be  no  less  than  (0.25//|  e | ) A. 

For  large  e,  the  restriction  on  the  lengths  of  the  antenna  can  be 
relaxed  so  long  as  they  are  situated  at  least  (0.7//[  e | ) X away  from 
the  ground  plane. 
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4.  ARBITRARY  ORIENTED  WIRE  ANTENNAS  OVER  IMPERFECT  GROUND 

In  this  section  expressions  for  mutual  impedance  between  two 
arbitrary  oriented  wire  antennas  over  plane  imperfect  ground  are  de- 
rived both  by  the  Sommerfeld  formulation  and  by  the  reflection- 
coefficient  method.  The  problem  configuration  consists  of  two  arbi- 
trarily oriented  current  elements  having  a current  distribution  I 

and  I . They  are  situated  over  a plane  imperfect  ground  as  shown 
b 

in  Fig.  10.  The  ground  is  characterized  by  the  complex  dielectric 
constant  e.  It  is  possible  for  1^  and  1^  to  overlap  one  another, 
thereby  forming  a wire  junction.  The  method  of  analysis  is  the  same 
as  presented  in  the  last  two  sections. 


4.1.  Exact  Sommerfeld  Formulation 


In  order  to  compute  the  mutual  impedance  between  two  arbi- 
trarily oriented  current  elements,  the  source  current  element  1^  is 

x y z 

split  up  into  three  components  1^,  1^  and  1^.  The  total  electric 

field  £ is  then  the  summation  of  ^ ^ and  E - fields  due  to 

IX  Iy  IZ 

A A A 

x y z 

the  component  current  elements  1^,  1^  and  1^  as  in  Fig.  10. 

Hence  from  (2.5),  (3.11),  (3.1A)  and  (3.15)  the  electric  field  due  to 

I,  can  be  written  as 
A 


[- 


Juhj 

A 7T 


— v~.  ][$  x x(« 

j (jJ A 7T c Jl  1A'80 


81  + 8sH^ 


+ 2 (82)ldtA 


(A 


Similarly  for 


Fig.  10.  Arbitrary  oriented  elements  over  an  imperfect 
ground  plane. 


-54- 


(4.2) 
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4ir 


°+3=!r1  » >*  (8o 


81  + gsH^ 


A l, 


+ 2 XA  9y  (g2^dilA 


and  from  (2.5),  (2.19)  and  (2.24)  the  field  is 


*,  • 1 
a a . 


% ■££-]  (2  l‘  (gQ  - 8,  + G„„)]dt, 


jw4ire 


sV' 


(4.3) 


Next  some  relations  are  presented  which  are  very  helpful  in  deriving  a 


formula  for  mutual  impedance  between  1 and  First 


3 f ,.33  , . 3 ,GSV, 

3x  (gsH)  3z  3x  g2  3x  ( e * 


3x'  K e } 
A 


3x'  (g3) 
A 


is  obtained  from  (3.22).  Also, 


3 r \j_9  3,,_3  /GsVx 

3y  + 3z  3y  (g2)  3y  ( E ) 


3 GsV,  9 , , 

K * 9yl  (83>  ; 


It  is  well-known  that 


3p 


(g0} 


- r-r  (gn)  » where  p may  be  'x',  ’y'  °r  ' z'  , 

3pA  0 


9 3 

(g,)  = - r-r  (g.)  , where  q may  be  V or  'y'  , 


3q  V61 

3 


K 


3q  ^gsH^  3q'  ’ 9q  ^2^  — ’ 


*\  _ I ' rtq  i 

9qA  2 


■h  <8i>  - si;  <*i>  “d  h <Gsv>  - sly  <V  • 

A A 


Applying  the  above  identities  the  total  E-field  can  be  obtained  in  the 
following  form  after  some  algebraic  manipulations  on  (4.1),  (4.2)  and 
(4.3). 
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uWta 
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$ 'I  V 


jojp 

4tt  ^A(g0  _ 81  + gsH^  + 2 {_IAQ  ” 


XA  * W}ld*A 


+ jco4iTG  IA  ‘ ^A(g0  81  + g3^  " TAT^d£/ 


(4.4) 


where  Q and  T are  given  by  the  following: 


[GsV  + 3z'  (g2^  _ 88H] 
A 


■ 2 Jn(Xp)  exp[-  V X 


— 

■-k,  (z  + z')]  i 


-V^>2 


(4.5) 


k2  v^2-k2  + k2  n2- k2 


T = 3z^  f2gl  GsV  g35  = g2  * kl 
A 


(4.6) 


Transferring  the  derivative  (^)  operating  on  the  Green's  function 
over  to  the  current  (I  ) and  assuming  the  current  to  be  zero  at  the 

A 

ends  of  the  open  wires, 


^A(g0  " 81  + gsH^  + 2 * XAQ  + d£  82^d£; 

A 


+ j oj 4 it e ^d£  (g0  ” 81  + g3)  ■ XAT^d^A 

o A 

A 

The  mutual  impedance  between  the  two  current  elements  "f  and  ^ 


(4.7) 


can  then  be  expressed  as 


= - t • f d£_ 

J A B B 


(4.8) 


After  some  algebraic  manipulations,  the  mutual  Impedance  reduces  to 


'AB 


4tt 
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dl, 


dV*A  * 81  + gsH)  + TB  {dl  82  " IAQ^ 
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J io4ttgo 


d*B 


dIA 

*di.  (g0  “ 81  + 83J 
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TZ  u2l 

XA  82kl 


(4.9) 


For  p < 0.5A,  the  semi-infinite  integrals  in  gg^»  g£,  and  Q are  com- 
puted in  an  essentially  exact  manner  along  a contour  shown  in  Fig.  2. 
For  p _>  0.5A,  a modified  method  of  steepest  descent  is  applied  to  take 
into  account  the  pole  near  the  saddle  point,  and  the  Integrals  can  be 
expressed  as  follows: 


Q = - g^  (cos  0 - \je  - sin^O)^ 
g2  * jg3  (cos  0 - \/e  - sin20) /k^ 


(A. 10) 
(4.11) 


and  ggH  and  g3  are  given  by  (3.29)  and  (3.30). 

A user-oriented  computer  program  has  been  written  and  is  presented  in  a 
later  report.  The  underlying  mathematics  of  the  program  is  relatively 
straightforward  once  the  general  formula  for  Z in  (4.9)  is  available. 

AJ5 


4.2.  Reflection-Coefficient  Method 

In  order  to  find  the  mutual  impedance  between  two  arbitrarily 
oriented  current  elements  over  the  plane  surface  of  an  imperfectly 
conducting  earth,  the  effects  of  the  direct  ray  and  reflected  ray  are 
considered  separately.  The  effect  of  the  direct  ray  is  straightforward 
and  is  the  same  as  in  [15].  In  order  to  obtain  the  effects  of  the  reflec- 
ted ray  a local  co-ordinate  system  centered  at  the  middle  point  of  dipole  A 
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is  chosen.  The  x and  z axes  of  the  co-ordinate  system  are  oriented 
In  such  a way  that  the  x-z  plane  passes  through  the  center-point  of 
dipole  B.  Also  the  x-y  plane  is  parallel  to  the  imperfectly  conduct- 
ing earth.  The  currents  of  element  A are  now  split  up  into  three 
components  each  being  parallel  to  the  local  co-ordinate  axes  chosen. 

As  seen  from  (2.49),  (2.50),  (3.40)-(3.42) , the  effect  of  the  reflected 
ray  is  equivalent  to  taking  the  original  current  element  and  placing  it 
at  the  image  point  C as  in  Fig.  10.  Then  the  contribution  to  the  reflected 

ray  due  to  I*  is  -I'TM(EX  + EZ  ) + r Ey  where  EX  , Ey  and  EZ  are  the 
a m _x  -X  l e,  _x  jX  jX  jX 

A A A A A A 

x 

three  components  of  the  electric  field  generated  by  1^.  Similarly  for 

Iy  the  corresponding  contribution  is  ^-^(E51  + EZ  ) + T _,Ey  and  for 

A r TM  y y TE  y 

A Lk  A 

Iz,  rrIM(EX  + EZ  ) + r„,_  Ey  . Hence  the  total  E-field  due  to  the 
A TM  z yz  TE  _z 

A A A 

reflected  ray  is  given  by 


r_,(*EX  + 2EZ  - St EX  - 2EZ  - StEX  - 2EZ  ) 

™ lz  lz  lx  lx  Iy  Iy 

A Lk  L\  Lk  A A 
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The  mutual  impedance  between  the  two  current  elements  I and  1^  is 


defined  as 
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where  Rx  - ^ ; R « yB  - and  Rz  = zB  + and  R2  =VRX  + Ry  + Kl 


The  total  mutual  impedance  between  the  two  current  elements  becomes 


Z = Z°  + Z1,  where  ZD  is  given  by. 
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where  is  given  by  (2.17). 

A user-oriented  computer  program  has  been  written  [17,24]  to 
solve  for  the  current  distribution  and  the  field  pattern  of  arbitrarily 
oriented  wire  antennas  (this  may  include  wire-junctions)  over  plane 
imperfect  ground.  Again  the  program  follows  directly  from  (4.13)  and 
(4.14). 

4.3.  Validity  of  the  Reflection-Coefficient  Method  as  Compared 

with  the  Exact  Solution 

In  general  it  is  difficult  to  predict  the  accuracy  of  the  so- 
called  exact  solution  below  a height  of  0.03X  from  the  ground  plane. 

If  sections  of  the  antenna  structure  below  this  height  are  all  vertical 
then  there  are  no  difficulties  but  if  there  is  any  horizontal  structure 
below  a height  0.03A  then  the  method  of  integration  used  to  integrate 
the  semi-infinite  integrals  fails.  It  is  safe  to  say  that  as  long  as 
there  are  no  horizontal  wire  elements  below  a height  of  0.03X  and  as 
long  as  any  structure  below  this  height  consists  of  all  parallel 
vertical  wires  then  the  exact  solution  is  valid. 

For  the  reflection-coefficient  method  it  is  impossible  to  make 
any  comment  regarding  accuracy  for  arbitrary  structure.  However,  as 
long  as  the  antenna  structure  is  within  the  limits  presented  in  the 
previous  sections,  it  is  possible  to  draw  conclusions  regarding  its 
relative  accuracy. 


[24]  T.  K.  Sarkar,  "Analysis  cf  Radiation  by  Arrays  of  Arbitrarily 
Oriented  Wire  Antennas  Over  Plane  Imperfect  Ground  (Reflection 
Coefficient  Method),"  IEEE  Trans.  Ant,  and  Propagat.  (to  be 
published) . 
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4.4.  Conclusion 


Two  user-oriented  programs  have  been  developed  using  both  the 
Sommerfeld  formulation  and  the  reflection-coefficient  method  based  on 
the  analysis  presented  in  the  last  section.  The  methods  presented  here 
and  their  corresponding  programs  handle  arbitrary  wire  geometry  includ- 
ing wire  junctions  and  use  piecewise  linear  functions  to  represent  both 

■a  and  V 


5.  OPTIMIZATION  METHODS  FOR  ARBITRARILY  ORIENTED  ARRAYS  OF  ANTENNAS 


i 


| 

i 

i 


f 


OVER  IMPERFECT  GROUND 


In  this  chapter  methods  are  presented  for  optimizing  certain  per- 
formance indices  of  arbitrarily  oriented  arrays  of  wire  antennas  over 
imperfectly  conducting  ground.  The  indices  considered  include  direc- 
tivity, maximum  power  gain,  quality  factor,  and  the  efficiency  index 
sometimes  referred  as  main-beam  radiation  efficiency.  Optimization 
problems  both  with  and  without  constraints  on  the  resulting  antenna 
pattern  can  be  handled.  In  the  first  six  sections  attention  is  re- 
stricted to  determine  array  feed  voltages  that  will  optimize  some  per- 
formance index.  Strait  and  Kuo  [25]  presented  the  optimization  methods 
for  arrays  of  parallel  vertical  wire  antennas  in  air.  An  attempt  has 
been  made  to  extend  these  methods  for  arbitrarily  oriented  structures 
now  situated  over  imperfect  ground.  The  next  section  deals  with  the 
optimization  of  certain  performance  indices  as  mentioned  but  this  time 
with  respect  to  the  geometry  of  the  antenna  structure.  This  includes 
optimization  with  respect  to  spacing  between  elements,  lengths  of  the 
elements  and  both  lengths  and  spacings.  Finally  optimization  of  the 
indices  with  respect  to  loads  applied  at  certain  points  on  the  antenna 
structure  is  also  considered. 

The  problem  of  interest  consists  of  an  arbitrarily  oriented 
antenna  array  of  N-input  ports  plus  a distant  test  antenna  having  one 
input  port  all  situated  over  an  imperfect  ground  plane.  This  is 
illustrated  in  Fig.  11.  The  array  and  the  test  antenna  form  an 


[25]  B.  J.  Strait  and  D.  C.  Kuo,  "Optimization  Methods  for  Arrays  of 
Parallel  Wire  Antennas,"  Scientific  Report  No.  18,  Contract  No. 
F19628-68-C-0180,  AFCRL-72-0725,  Syracuse  University,  Syracuse, 
New  York:  December  1972. 
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(N+l)-port  network  whose  terminal  characteristics  are  determined  by  an 
(N+l)  by  (N+l)  matrix.  In  order  to  obtain  the  terminal  characteristics 
the  antenna  array  is  divided  into  subsections  (as  in  the  "Moment  Method") 
and  the  current  of  the  test  antenna  is  represented  by  one  expansion  func- 
tion only.  The  generalized  impedance  matrix  is  computed  which  completely 
characterizes  the  wire  structure.  Either  the  "exact"  solution  or  the 
reflection-coefficient  method  is  used  to  calculate  individual  matrix  ele- 
ments but  one  method  is  adhered  to  compute  the  whole  matrix.  The  formulas 
presented  in  Chapters  2,  3 and  4 are  utilized  depending  on  the  method  used 
and  the  nature  of  the  antenna  structure.  The  impedance  matrix  is  inverted 

to  obtain  the  generalized  admittance  matrix  f Y ] of  the  entire  structure. 

T 

The  terminal  admittance  matrix  [Y]  is  obtained  by  retaining  only  those 
elements  Y for  which  both  i,j  represent  the  input  ports  of  the  antenna 


structure.  So, 


I 

r 

t 

[Ia] 
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l\al 

tt 
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V 

t 

jv 

(5.1) 


where  and  Y are  the  terminal  current,  voltage  and  the  input  admit- 

tance of  the  test  antenna,  respectively.  The  corresponding  terminal  char- 
acteristics for  the  antenna  array  are  expressed  by  the  matrices 
[V  ]„  , and  [Y  ]„  [Y  „ and  [Y  ]..  , represent  the  mutual  admittances 
between  the  test  antenna  and  the  array  ports  incorporating  the  effects  of 

^ /V 

the  imperfect  ground  plane.  Also  [Y  ] = [Y  ] and  [Y  ] = [Y  ] . Once  the 

aa  aa  ta  at 

four  admittance  matrices  of  (5.1)  are  known  array  optimization  problems  can 
be  handled  as  illustrated  in  the  next  sections. 

5.1.  Gain  Maximization 


The  power  gain  for  an  antenna  array  is  defined  by 
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^ _ Att  x (Radiation  intensity  for  the  specified  direction) 
Power  input  to  the  array 


(5.2) 


The  power  input  to  the  array  is  given  by 


Pln  - j [V  *1  {[Y  ) + [Y*  ]}  [V  J 

in  2 a aa  aa  a 


(5.3) 


where  * denotes  complex  conjugate  and  ^ denotes  the  transpose  of  the 
matrix.  The  test  antenna  is  assumed  to  be  appropriately  polarized  and 
located  in  the  specific  direction  for  which  G is  desired.  It  is  also 
assumed  that  the  test  antenna  is  at  a distance  r from  the  array  and 
situated  in  the  far  field.  Under  polarization-matched  conditions,  the 
square  of  the  magnitude  of  the  test  antenna  terminal  current  is  propor- 
tional to  the  radiation  intensity  (R)  of  the  incident  field.  For  V = 0, 
the  power  received  (P  ) is  related  to  the  gain  (g)  and  the  receiving 
aperture  (A)  of  the  test  antenna  by 


111-  R. 


tt 


(5.4) 


where  R is  the  input  impedance  of  the  test  antenna.  Since  g,  A and 
' - the  wavelength  of  propagation,  are  related  by 


A = AjL 
4 71 


it  follows  that 


2 r2  iZ  2 K 

*-  -T7—  "t1  & 

4tt  A 

K is  a constant  which  depends  solely  on  the  test  antenna.  Since 

I = [Y  ][V  ] for  V = 0,  (5.6)  can  be  expressed  as 
t ta  a t r 

R - j-  [V*]  [7*  ] [Y  ] [VI 

ott  a ta  ta  a 


(5.5) 


(5.6) 


(5.7) 
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Application  of  (5.3)  and  (5.7)  to  (5.2)  yields 


G 


Iv;>  [Y?.l  <V 


lV.*l  <tVM]  + [Y*,])  [Va) 


(5.8) 


The  procedure  for  optimizing  (5.8)  for  some  specified  direction  is 

well-known  [26].  This  is  a ratio  of  Hermitian  quadratic  forms,  and  it 

follows  from  the  properties  of  these  network  matrices  that  { [Y  ] + [Y*  ]} 

aa  aa 

is  positive  definite  while  [Y*  ] [Y^  ] is  positive  semidefinite.  Hence 

ta  ta 

all  eigenvalues  associated  with  the  following  equation 


|Y?a]  'V  >V-Ft|VJ  + [Yaa])  'V 


(5.9) 


are  either  zero  or  real  and  positive.  Since  [Y*  ] [Y  ] is  a one-term 

v ta  ta 

dyad  there  is  only  one  non-zero  eigenvalue,  and  the  corresponding  eigen- 
vector is  found  from 


- mi  + [T*  j} 

G aa  aa 
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[Y?aJ  |Ytal  '''a1 


'V 


(5.10) 


Since  [Y  ] [V  ] is  a scalar  the  excitation  voltages  required  to  maximize 
ts  a 

the  gain  in  the  specified  direction  are  contained  in  the  eigenvector 


[V  ] = C {[Y  ] + [Y*  ]}_1  [Y*  ] 

a aa  aa  ta 


(5.11) 


where  C is  a constant.  Substitution  of  (5.11)  to  (5.10)  yields 
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max 


K[Y  ] 
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ir1  [y*  i 

ta 


(5.12) 


[26]  R.  F.  Harrington,  "Field  Computation  by  Moment  Methods,"  The 
Macmillan  Company,  New  York,  1968. 
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Hence  the  voltages  given  by  (3.11)  (the  constant  C can  be  dropped  with 

no  effects  on  the  result)  can  be  used  to  maximize  the  gain  of  the  array 

in  the  given  direction.  The  result  is  G given  by  (5.12). 

° max  ° 


5.2.  Efficiency  Indices 


An  efficiency  index  S (or  inverse  of  sensitivity  factor)  sometimes 
called  the  main-beam  radiation  efficiency  [27]  has  been  defined  for  arrays 
of  isotropic  point  sources  as 

g _ (radiation  intensity  corresponding  to  the  direction  of  max,  radiation) 

(sum  of  the  excitation  current  magnitudes  squared) 

(5.13) 

to  serve  as  a measure  or  indication  of  the  supergain  condition.  It  has  been 

pointed  out  [27]  that  superdirective  arrays  may  require  very  large  currents 

of  opposite  signs  in  neighboring  elements,  resulting  in  excessive  heat  loss 

and  very  low  radiation  intensity  in  the  direction  of  the  main  beam.  This 

generally  undesirable  condition  is  indicated  by  a relatively  low  value  of  S. 

The  efficiency  index  has  similar  significance  for  arrays  of  electrically 

short  wires.  Use  of  (5.7)  and  [I  ] = [Y  ] [V  ] in  (5.13)  yields 

a aa  a 

K IV*]  [Y*  1 [Y  1 [V  ] 

S = £— - — — — (5.14) 

IV*]  [Y*  ] [Y  ] [V  ] 
a aa  aa  a 

For  longer  wires  it  may  be  more  useful  to  deal  with  either  one  of  two 
possible  alternative  quantities  that  will  also  be  labeled  efficiency 
indices.  One  of  these  is  obtained  by  using  (5.13)  with  its  denominator 
replaced  by  the  sum  of  the  excitation  voltages  squared.  By  denoting  the 
resulting  quantity  by  and  use  of  (5.7)  in  (5.13)  yields 


[27]  D.  K.  Cheng,  "Optimization  Techniques  for  Antenna  Arrays,"  Proc. 
IEEE,  No.  12,  pp.  1664-1674,  December  1971. 
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(5.15) 


v [V*]  [Y*  ) [Y  ] [Vl 

c — iS  ^ tfl  t3  d 

1 [V*]  [V  ] 

a a 


For  longer  wires  a suitable  alternative  is  obtained  by  using  (5.13)  with 
its  denominator  replaced  by  the  sum  of  the  current  magnitudes  squared 
over  the  entire  array.  The  new  efficiency  index  is  then  obtained  as, 


s K 1V:»  |YL]  |Yta' 

2 2 [V*]  |Y*]  [Y]  [V  ] 

a a 


(5.16) 


where  the  admittance  matrix  [Y]  is  obtained  by  retaining  only  those 
columns  of  [Y]  which  correspond  to  the  input  ports  of  the  antenna 
array. 

The  excitation  voltages  that  will  maximize  the  efficiency  indices 
S,  and  can  be  obtained  from  the  method  of  Section  5.1.  Results  are 


S = | [Y  ] f(Y*  ] [Y  l)*1  [Y*  ] 
max  2 ta  aa  aa  ta 


(5.17) 


corresponding  to  the  feed  voltages  given  by  (within  an  arbitrary  constant) 


|Va]  ' <^a>  'V''1  I7?.) 


(5.18) 


Also , 


(S,)  = ^ [Y  ] [Y*  ] 

1 max  2 ta  ta 


(5.19) 


with  feed  voltages  given  by 


'V  ‘ ,Yta> 


(5.20) 


Finally 


<VMx  ■ f 'ha1  "Y*]  mr1  IV*aJ 


resulting  from  feed  voltages  given  by 


(5.21) 
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(5.22) 
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Thus  the  three  choices  for  efficiency  index  can  be  calculated  for  any 
arbitrarily  oriented  antenna  array  over  an  imperfect  ground  plane  by 
using  the  formulas  for  impedances  described  in  Chapters  2,  3 and  A. 


5,3.  Quality  Factor 

A quality  factor  can  be  defined  for  arbitrarily  oriented  arrays 
over  imperfect  ground  that  relates  gain  and  efficiency  index.  For  arrays 
of  electrically  short  wires  it  is  convenient  to  use  the  quality  factor  Q 
defined  for  arrays  of  isotropic  point  sources  as 


Ait  (sum  of  excitation  currents  squared) 
^ power  input  to  the  array 


(5.23) 


So  that  from  (5.2),  (5.13)  and  (5.23) 


G • QS 


(5.2A) 


Application  of  (5.3)  and  [I  ] = [Y  ] [V  ] in  (5.23)  yields 

a S3  3 

8n[V*  ] [Y*  ] [Y  ] [V  ] 

„ a aa  aa  a 


[V*J  UY  1 + [Y*  )}  [V  ] 
a aa  aa  a 


(5.25) 


Similarly,  alternative  quality  factors  and  are  defined  correspond- 
ing to  the  efficiency  indices  of  (5.15)  and  (5.16).  The  result  is  that 
(5.2A)  holds  with  the  product  QS  replaced  by  Q^S^  or  Q2S2»  which  ever  is 
appropriate.  By  using  the  equivalent  of  (5.2A),  then  it  is  obvious  that 


^1  = 
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(5.26) 
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and 


8ti  (V*  ] [ Y*  1 [Y]  IV 

a a 


(5.27) 


[V* | f I Y I + [Y  * 1 } [V  ] 
a aa  a a a 


Hence  the  quality  factors  <>,  <)  and  can  be  computed  once  the  antenna 
structure  and  the  feed  voltages  are  known.  These  quantities  can  be  opti- 
mized bv  a procedure  as  outlined  in  the  previous  section. 


5,4.  Opt imizat  ion  Subject  to  Oonst raints  on  Pattern  Nulls  in  the 
I'pper  Half-Space 

The  problem  of  optimizing  t he  directivity  or  any  of  the  efficiency 
indices  subject  to  constraints  on  the  resulting  pattern  nulls  in  the  upper 
half-space  is  discussed  in  this  suction.  The  field  pattern  in  a certain 
direction  above  the  ground  plane  for  an  arbitrarily  oriented  antenna  array 
over  imperfect  ground  will  be  proportional  to  | , where  f is  the  short- 

circuit  current  flowing  through  the  test  dipole  located  in  that  particular 
direction.  Specifically  the  field  (i’.r ) in  direction  p above  ground  is 
related  to  the  short-circuit  current  in  the  test  dipole  which  is  located 
iii  the  p-direction  by 

i K*’!  = /2nRtt  i ij’l  (5.28) 

where  ” and  K are  the  wave  impedance  and  input  resistance  of  the  test 
dipole  respectively.  The  problem  is  to  determine  the  feed  voltages  that 
will  provide  pattern  nulls  in  p specified  directions  in  the  upper  half- 
space and  a maximum  of  one  of  the  performance  indices  (say  the  gain)  in 
a given  (p+1)  direction  subject  to  these  constraints.  An  equation  re- 
lating the  short  circuit  current  in  the  test  dipole  for  each  of  the  (p+1) 
directions  to  the  admittance  matrix  and  the  feed  voltages  can  he  written 
as  follows: 
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where  the  matrix  [Y^]  is  the  transfer  admittance  matrix  between  the  test 
dipole  in  the  p-direction  and  the  excited  ports  of  the  antenna  array, 
taking  into  account  the  effects  of  the  imperfect  ground  plane.  i^~nR  " |l^| 
is  then  the  field  in  the  direction  designated  for  the  maximum  performance 
index.  Equation  (5.29)  can  be  portioned  in  the  following  form 


Hence 


[0] 

1 

s 

L * J 

Iva' 

ss  — 

[Y11] 
t Jp*P 


[Y21] 

1 t Jl*p 


rj1]-1  [Yt12] 


[Y12] 

1 t Jpx (N-p) 


[Y22] 

Ut  Jlx(N-p) 


l2}  A [B][V2] 
a — a 


tVa]pxl 


[V(N-p)xl 


(5.30) 


(5.31) 


and 


It  = ([Y21]  [B]  + [Y22]}  [v2]  A [C]  [v2] 


(5.32) 


If  the  performance  index  chosen  to  be  maximized  is  the  gain  then  from 
(5.2),  (5.6)  and  (5.32) 


G = 


K [v2*3  [c*]  [c]  iv2] 

^ a o 

2 P. 

in 


(5.33) 


and  P.  from  (5.3)  is  given  as 
In 
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ft /IV, to  1 >Vr,  ^ 


^ w » w ■<  h»h t*ftur»+  ww 


" ’SPfM*to».v-i«K«v  !•■•»*.  ■■ 


Pl„  ‘ 2 (VS]  fl\J  + [Vaa'  > 'Val 


and  with  [V  ] written  in  partitioned  form 

d 


Pi„-7t[Va‘):'Va*1,{1Yaa1+[^a1  1 


1^1 


[V«] 

& 


(5.34) 


Application  of  (5.31)  leads  to 


Pin-I  1Va‘l  UB*1  ! [U]>  {[Yaa]  ♦ [Ya*  ] } 


or 


P,n  iy  [v2*]  [Q]  [v2] 

in  L a a 


where  [U]  is  a unit  matrix  and  [Q]  is  defined  by  (5.34)  and  (5.35). 
Hence,  application  of  (5.35)  to  (5.33)  yields 


(5.35) 


G = K 


[V?*]  [7*]  [C]  [V2] 

a a 

~2*  2 

[Vj  ] [Q]  [Vf] 

o a 


(5.36) 


Once  again  (5.36)  is  a ratio  of  Hermitian  forms  with  [Q]  positive  definite 
and  [C*]  [C]  positive  semidefinite.  Furthermore  the  latter  is  a one  term 
dyad  so  that  the  procedure  for  maximizing  the  gain  is  identical  to  that 
used  in  Section  5.1.  The  results  are 


G * K • [C]  [Q]*1  [C*] 
max 


(5.37) 


and 


[V2J  = [Q]"1  [?] 


(5.38) 

2 1 

Once  [V  ] is  known  [Vi]  can  be  obtained  from  (5.31)  which  completes  the  solu- 
a a 

tion.  If  no  pattern  constraints  are  used  then  (5.37)  and  (5.38)  reduce  to 
the  solution  of  Section  5.1. 
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The  procedure  for  optimizing  any  of  the  efficiency  indices  defined  in 


Section  5.2  subject  to  constraints  on  the  resulting  pattern  nulls  in  the 
upper  half-space  follows  exactly  the  same  procedure. 

5.5.  Optimization  Subject  to  Constraints  on  the  Resulting  Sidelobes  in 
the  Upper  Half  Space 


Next,  the  problem  of  finding  excitation  voltages  that  will  maximize 
a performance  index  (of  an  arbitrarily  oriented  antenna  structure  over  an 
imperfect  ground  plane)  in  a given  direction  subject  to  constraints  on  the 
resulting  sideline  levels  is  considered.  Suppose  1°  denotes  the  short 
circuit  current  in  a small  test  antenna  located  in  a direction  corresponding 
to  the  direction  designated  for  a maximum  of  a given  performance  index  (say 
gain),  and  1^,  1^,  ...»  Ip  denote  values  of  the  short  circuit  current, 
induced  in  small  fast  antennas  located  in  direction,  corresponding  to  the 
peaks  of  the  p sidelobes  to  be  constrained  in  the  pattern.  The  short-circuit 
currents  in  the  test  antennas  are  proportional  to  the  field  intensities  in 
those  directions  and  are  related  by  (5.28).  Using  partitioned  matrices 


r 

1 1° 

_t_ 

lTt  1* (N-p) 

[yAB] 

1 t Jl*p 

x1 

t 

« 

• 

= 

— - - - - 

- - - 

[VaJ(N-p)xl 

|VaW 

Ip 

t 

_ 

[y®a] 

t px(N-p) 

[yrh] 

t Jp*p 

BB 

Assuming  [Y  ] to  be  nonsingular 
t "Kl-r 


[V®]  = [Y®8]'1 


|-  [y®a]  [VA] 
t a 


- 


(5.39) 


(5. AO) 
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An  alternative  to  (5.40)  is  possible  if  it  is  not  required  to  satisfy  (5.39) 

exactly.  That  is,  if  a least  squares  solution  is  acceptable  it  is  possible 

B BA 

to  use  less  than  p voltages  in  [V  ].  That  is,  [V  1 and  (V  I are  of  dimen- 

a a a 

sions  m and  N-m  respectively  where  m < p,  then  the  dimensions  of  the  other 

BB 

matrices  change  accordingly  with  [Y^  ] no  longer  square.  A least  squares 
solution  is  then  obtained  using 

[V®]  - {[YBB*]  [Y88]}"1  [Y8B*] 

a Cl  L 

in  place  of  (5.40) 


r£] 


[yba]  [va] 


(5.41) 


Using  a suggestion  of  Sanzgiri  and  Butler  [28]  the  short-circuit 
current  in  the  test  dipole  corresponding  to  the  peak  of  the  ith  sidelobe 
can  be  denoted  by  I*  = 1°  so  that 


Then,  from  (5.38) 


where 


and 


ft1" 

,0 

e , I 

t 

1 t 

• 

• 

• 

a 

• 

• 

Ip 

e 1° 

t 

P t 

* = 1 - [Y^]  [YBB]_1  [E] 

[F]  = [Y^]  - [YAB]  [YBB]_1[YBA] 


(5.42) 


(5.43) 

(5.44) 


(5.45) 


[28]  S.  M.  Sanzgiri  and  J.  K.  Butler,  "Constrained  Optimization  of  the 

Performance  Indices  of  Arbitrary  Array  Antennas,"  IFEE  Transactions 
on  Antennas  and  Propagation,  Vol.  AP-19,  No.  4,  pp.  493-498, 

July  1971. 


-74- 


The  input  power  is  given  by  (5.34)  and  using  (5.40)  and  (5.42)  - (5.45) 
in  (5.34)  results  in 


p.  = 7 ivf  ] [Q*]  [/] 

in  2 a a 


(5.46) 


where  in  this  case 


[Q * ] = {[U]  ; [G*]}  { [Y  ] + [Y*  ]} 
t aa  aa 


[G]  = ^ [Y®8]-1  {[E]  IF)  - 0[Y8A) } 


(5.47) 


(5.48) 


Of  course  (5.44),  (5.45),  (5.47)  and  (5.48)  must  be  changed  appropriately, 
if  (5.41)  is  used  rather  than  (5.40). 

If  the  performance  index  to  be  maximized  is  chosen  again  to  be  gain 


then , 


G - K 


[Vfl  [F*]  [F]  [V^] 

ft*]  IQ’] 


(5.49) 


The  matrix  [ F* ] [F]  is  a one  term  dyad  and  the  procedure  for  optimizing 


(5.49)  is  again  that  of  Section  5.1.  The  results  are 

G„a*  ■ rW  miQ’fhr.) 


[/]  = iqT1  [F^] 

a 


[V8)  = [G]  [/] 


(5.50) 


(5.51) 


(5.52) 


The  iterative  optimization  procedure  is  then  simply  to  guess  the 

initial  sidelobe  directions  and  compute  the  matrices  necessary  to  determine 

G and  the  required  voltages  (5.50)  to  (5.54).  The  resulting  sidelobe 
max 


directions  and  peaks  are  then  computed  and  compared  with  desired  results. 
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The  new  sicilobe  directions  are  used,  of  course,  at.  "initial  guesses" 
for  the  next  iteration.  The  iterative  procedure  is  continued  until  the 
levels  of  all  sidelobes  to  be  constrained  are  within  some  prescribed 
tolerance. 

Obviously,  constrained  optimization  of  the  various  erficiency  Indices 
defined  earlier  can  be  carried  out  in  the  same  way. 


5.6.  Optimization  of  One  Performance  Index  Subject  to  a Constraint 
on  Another 

A procedure  for  optimizing  a given  performance  index  (of  an  arbi- 
trarily oriented  antenna  array  over  an  imperfect  ground  plane)  subject  to 
a constraint  on  another  performance  index  is  outlined  here.  In  this  dis- 
cussion an  index  denoted  by 


"v 

"V  <v 

will  be  optimized  subject  to  the  constraint  given  by 


■fo  "V  lv,' 

1^1  IV  <V 


= Y (a  given  constant). 


(5.53) 


Once  the  performance  indices  are  specified  [ R ^ ] , [J^l.  an<*  have 

precise  meanings.  A solution  is  obtained  by  setting 


, |v;'  ‘V  |v,]  , . 
in1  "v  'V 


1 [V*]  [R,l  IV  ] 

a J a 

[V*]  [R  ) [V  ] 
a 4 a 


- Y 


(5.54) 


stationary  with  respect  to  cne  column  vector  [V  ] and  the  Lagrangian 


multiplier  >.  Setting  the  first  variation  of  L with  respect  to  [V  ] 

3 

equal  to  zero  yields 
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(5.55) 


[V  1 = a IH]"1  [S*] 
a i 

where  a is  a constant  given  by 


{ [sxl 


[V  ]} 
a 


Uv*] 

a 


[R21 


[V  ]) 
a 


{[V*] 

a 


[Rj] 


[V]} 

a 


-1 


and 


[1^]  = (s*l  ls1) 


(5.56) 


The  matrix  Ini'!  gi'en  by 


IH]  = IR2]  + yblRj  - b[R3]  (5.57) 

where  the  constant  b is 

b = > ([V*]  [R  ] [V  ]}_1  { IV*)  [R  ] [V  ]}2  {[V*]  [Rn  ] [V  J)'1  (5.58) 

a 4 a a 2 a ala 

The  results  are  essentially  those  obtained  by  Lo,  Lee  and  Lee  [29]  for  the 
special  case  when  [R^]  is  a unit  matrix;  although  of  course,  not  only  mutual 
coupling  is  incorporated  in  the  matrices  dealt  with  here,  but  also  the  effects 
of  the  imperfect  ground. 

When  (5.55)  is  substituted  in  the  constraint  equation  (5.53)  it  follows 

that 


[Sj  UH]"1  (y(R41  - [R3D  [H]"1  [S*]}=  0 (5.59) 

where  it  is  assumed  that  I ] > [ R ^ ] , [R^]  and  [H]  are  Hermitian.  The 
only  unknown  in  (5.59)  is  b which  is  contained  in  matrix  [H]. 

The  procedure  for  finding  b is  spelled  out  by  Winkler  and  Schwartz 
[30].  Tt  is  evident  from  (5.59)  that  the  column  vector  shown  in  brackets 


[29]  Y.  T.  Lo,  S.  W.  Lee  and  D.  H.  Lee,  "Optimization  of  Directivity  and 
Signal  to  Noise  Ratio  of  an  Arbitrary  Antenna  Array,"  Proc.  IEEE, 
Vol.  54,  No.  8,  pp.  1033-1045,  August  1966. 

[30]  L.  P.  Winkler  and  M.  Schwartz,  "A  Fast  Numerical  Method  for 
Determining  the  Optimum  SNR  of  an  Array  Subject  to  a Q-Factor 
Constraint,"  IEEE  Transaction  on  Antennas  and  Propagation,  Vol. 
AP-20,  No.  4,  pp.  503-505,  July  1972. 
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is  orthogonal  to  [S^].  It  follows  that  a complete  set  {[f^]}  can  be  con- 
structed with  [S1]  as  one  of  its  elements.  If  [S^]  is  written  as 

[S^]  = [s^,  (5.60) 

where  N is  the  number  of  independent  applied  excitations  then 

[fj]  — [Sj*  s^,  S ^ • t 

I f o ] = f r1  . 0,  o o] 

1 1 2 

Ifi]  = r- » °*  7"’ 0 0] 

. 3 


IfJ  = [-  0,  0,  0, 

S1 


(5.61) 


forms  a complete  independent  set  with  | f ^ ] orthogonal  to  all  others.  The 
bracketed  column  vector  in  (5.59)  can  then  be  expressed  as  a linear  combi- 
nation of  [f  ] written  as 

^ N ^ 

inf1  (t|K4J  - [ «3 1 ) IH]"1^*]  = l q1I?1J  (5.62) 


This  is  easily  arranged  as 

[w]  [o']  = 0 (5.63) 

where 

IQ  1 = (-  i»  q2»  q3.  *••*%]  (5.64) 

and  [W]  is  a matrix  with  columns  [ ] given  by 

[Wj]  = [7*]  = [S*J  (5.65) 

[w± ] = h2 (y(R4  1 - [R3D  \T{]  +?b[R2J[f1] 

+ [R21  (y(R41  - fR3])-1  [R2 ] [Tj] 

for  i = 2,  3,  4,...N  (5.66) 
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Rewriting  (5.63)  in  terms  of  its  real  and  imaginary  parts  in  the  form 


or 


— — 

— — 

w - w w - u 

Hr  lli  w12r  12i  ' * * 
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q2i 

• 

• 
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«■»  — . 

_ -J 

[W]  IQ']  = 0 


(5.67) 


(5.68) 


where  [W  ] is  a square  matrix  of  dimension  twice  that  of  [W],  and  [cf]  is 
a column  vector  of  dimension  twice  that  of  [Q] . By  using  (5.66),  (5.68) 
can  be  rewritten 


(b  [A^ ] + b [A2]  + [ ] } [Q ' ) = 0 


(5.69) 


where  [ A^  ] and  [A^  ] are  2N  * 2N  singular  matrices  (their  first  two  columns 
are  zero)  and  [A3]  is  a 2N  x 2N  invertible  matrix.  Then  from  (5.69) 

Winkler  and  Schwartz  [30]  show  that  values  of  b can  be  found  by  determining 
the  real  eigenvalues  of 


[G]  [X]  = i [X] 

where  [G]  is  an  unsymr  . ric  matrix  given  by 

[0] 


[G] 


[U] 


- [A3]'1  [Aj] 


- [A3]”1  [A2] 


4N*4N 


(5.70) 


(5.71) 


where  [U]  is  a unit  matrix. 
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For  each  real  value  of  b one  set  of  excitation  voltages  is  computed 
using  (5.55).  Each  of  these  sets  is  then  used  to  compute  corresponding 
values  of  both  the  constrained  index  and  the  index  to  be  optimized.  The 
set  finally  chosen  is  the  one  which  yields  the  optimum  result  subject  to 
the  correct  constraint. 

5.7.  Optimization  with  Respect  to  Antenna  Geometry 

With  an  appropriate  optimization  procedure,  the  performance  indices 
presented  in  the  previous  sections  may  be  maximized  through  adjustments 
of  the  antenna  structure.  In  this  case  either  lengths,  spacings  or  both 
lengths  and  spacings  are  to  be  determined  so  as  to  maximize  a performance 
index  given  by  either  (5.8),  (5.14),  (5.15),  (5.16),  (5.25),  (5.26)  or 
(5.27).  The  starting  point  in  the  optimization  procedure  can  be  an  initial 
state  corresponding  to  some  classical  or  other  improved  design  procedure. 

At  each  refinement,  the  admittance  matrices  [Y  ] and  [Y  ] are  recomputed 

ucl  33 

and  the  performance  index  is  calculated  by  the  optimization  procedure  till 
the  design  parameters  meet  certain  error  specifications. 

The  next  problem  is  to  reduce  the  sidelobe  levels  below  a certain 
level  of  an  antenna  array  over  an  imperfect  ground  plane  by  adjusting 
the  antenna  geometry.  Since  the  electric  field  in  the  upper  half  space 
in  a certain  direction  is  related  to  the  short  circuit  current  1^  in  the  test 
dipole  by  (5.28),  an  error  criterion  can  be  defined  as 


n 

i i 1 2 i 1 2 

for 

i i | 
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= l 

i-1 

htl  - htl 

i’ti  * 

= 0 

for 

idi  < 

i*,i 

t 
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where  I is  the  value  of  the  short  circuit  current  in  the  test  dipole 
oriented  in  the  i-th  direction  and  n is  the  number  of  the  observations 
calculated  over  the  sidelobe  region  of  the  pattern.  Of  course,  the  error 
criterion  e is  the  quantity  to  be  minimized  as  the  optimization  procedure 
progresses.  Since  at  each  optimization  step  the  admittance  matrices 
[Y  ] and  [Y  ] are  recomputed  this  is  indeed  a very  time-consuming 
process.  These  types  of  optimizations  have  been  done  for  parallel  verti- 
cal wires  in  echelon  [31],  and  could  easily  be  extended  to  arrays  of 
arbitrarily  oriented  wires  over  imperfect  ground  provided,  of  course,  that 
sufficient  computer  time  is  available. 

5.8.  Optimizations  of  Loaded  Antennas 

In  this  problem  (as  in  the  previous)  one  has  to  recompute  the 

admittance  matrices  [Y  ] and  [Y  ] at  each  optimization  step.  The 

t3  33 

problem  in  this  case  involves  optimization  of  the  loads  that  are  applied 
either  at  some  specific  points  or  distributed  along  a section  of  the 
arbitrarily  oriented  antenna  structure.  Once  the  total  impedance  matrix 
[ Z ] ^ for  the  entire  structure  has  been  computed  excluding  the  loads  (but 
including  the  effects  of  the  imperfect  ground  plane)  the  effect  of  the 
loads  can  be  taken  into  account  by  adding  a load  impedance  matrix  [Z]^  to 
[Z]s<  This  load  impedance  matrix  is  a square  matrix  having  the  same 
number  of  columns  as  [Z]g,  where  the  number  equals  the  number  of  segments 
into  which  the  N-port  antenna  array  has  been  divided.  For  passive  point- 
by-point  loading  [ Z ] ^ is  a diagonal  matrix  having  a zero  on  the  diagonal 

[31]  B.  J.  Strait,  T.  K.  Sarkar  and  D.  C.  Kuo,  "Special  Programs  for 
Analysis  of  Radiation  by  Wire  Antennas,"  on  Contract  F19628-73-C- 
0047,  AFCRL-73-0399 , Syracuse  University,  Syracuse,  NY  13210. 
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corresponding  to  each  unloaded  triangle  function.  Thus  if  there  is  a 
lumped  load  of  Z ohms  at  a wire  position  corresponding  to  the  peak  of 
the  nth  triangle  function,  then  the  element  corresponding  to  the  nth 
row  and  column  of  the  diagonal  matrix  [ Z ] ^ is  Z^.  Thus  the  total  im- 
pedance matrix  [Z]^  comprising  of  the  louded  antenna  structure  over  the 
plane  imperfect  ground  is  then 


1 


[Z]T  = [Z]s  + [Z]l 

T 

In  this  case  the  terminal  admittance  matrix  [Y]sf  is  obtained  by  retain- 
ing only  those  elements  {[Z]^  which  correspond  to  the  excitation 

ports  of  the  antenna  structure.  Hence  once  [Y  ] „ and  [Y  ] „ are  ob- 

ta  si  aa  s£ 

tained  the  performance  indices  given  by  either  (5.8),  (5.14),  (5.15), 
(5.16),  (5.25),  (5.26)  or  (5.27)  can  be  easily  computed.  An  optimization 
procedure  now  can  be  used  with  an  initial  guess  to  optimize  the  perfor- 
mance index  with  respect  to  the  loads  Z . The  starting  point  for  the 
iterative  procedure  is  the  unloaded  structure.  The  final  result  is  a 
structure  loaded  at  n-points  by  [ZJn]  which  optimize  the  required  per- 
formance index. 

The  same  procedure  can  again  be  applied  to  reduce  sidelobes  of 
the  fielos  in  the  upper  half  space  by  application  of  loading  to  antenna 
structure.  As  described  in  the  previous  section,  the  same  error  cri- 
terion is  defined  but  the  parameters  to  be  optimized  are  not  the  struc- 
ture over  the  plane  imperfect  ground.  The  starting  point  for  the  iterative 
procedure  is  the  pattern  in  the  upper  half  space  of  an  unloaded  antenna 
geometry.  The  final  result  is  the  same  array  over  the  plane  imperfect 
ground  but  loaded  and  having  sidelcbe'  near  the  desired  level  in  the  upper 
half  space. 
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6.  CONCLUSION 


The  mathematical  foundations  for  newly  developed  user-oriented 
computer  programs  have  been  presented  and  described  for  handling  arbi- 
trarily oriented  thin-wire  antennas  radiating  in  the  presence  of  a 
plane  imperfect  ground.  Properties  of  the  imperfect  ground  are  taken 
into  account  either  exactly  by  using  the  Sommerfeld  formulation  or 
approximately  by  using  the  reflection-coefficient  method.  Application 
of  interpolatory  quadrature  formulas  to  integrate  the  semi-infinite 
integrals  encountered  in  the  exact  Sommerfeld  formulation  has  reduced 
the  time  of  computation  as  compared  with  other  available  procedures. 

A modified  method  of  steepest  descent  has  been  used  to  evaluate  the 
semi-infinite  integrals  to  reduce  required  computing  time  without 
significant  loss  of  accuracy.  The  commonly  used  reflection-coefficient 
method  is  derived  in  detail  and  its  relative  accuracy  is  discussed. 

Methods  have  been  presented  fox  optimizing  various  performance 
indices  of  arbitrarily  oriented  thin-wire  antennas  over  plane  earth 
ground.  Performance  indices  considered  include  directivity,  maximum 
power  gain,  quality  factor  and  efficiency  index.  It  has  b?en  shown 
how  a performance  index  can  be  optimized  subject  either  to  constraints 
on  the  directions  of  nulls  in  the  resulting  pattern  or  to  constraints 
on  the  levels  of  the  resulting  sidelobes.  It  has  also  been  shown  how 
one  performance  index  can  be  optimized  subject  to  a constraint  on 
another.  The  methods  presented  are  quite  general  in  that  the  wires 
ran  be  excited  or  loaded  at  arbitrary  points  along  their  lengths. 

In  summary  the  contributions  presented  in  this  thesis  include 
the  fol lowing: 


-83- 


1)  Application  of  interpolatory  quadrature  formulas  to  Integrate  the 
semi-infinite  integrals  encountered  in  the  exact  Sommerfeld  formulation 

was  illustrated.  » 

2)  Use  of  a modified  method  of  steepest  descent  to  evaluate  the  semi- 

t 

Infinite  integrals  to  reduce  the  time  of  computation  without  significant 
loss  of  accuracy  has  been  demonstrated. 

3)  The  commonly  used  reflection-coefficient  method  was  derived  in  de- 
tail and  its  relative  accuracy  was  discussed, 

4)  New,  user-oriented  computer  programs  for  treating  the  imperfect 
ground  problem  have  been  developed  and  described. 

5)  Formulas  and  results  convenient  for  application  to  optimization 
and  design  problems  have  been  presented. 


ir 


Appendix  A.  QUADRATURES  OF  THE  HIGHEST  ALGEBRAIC  ACCURACY 

fh  r,n 

Certain  integrals  of  the  type  f(x)dx  and  F(x)dx  have  been 

a 0 

encountered  in  the  course  of  this  work.  A method  for  evaluating  these 
types  of  integrals  very  efficiently  is  discussed  here.  The  following 
subject  matter  has  been  summarized  from  Krylov  [33]  to  illustrate  the 
salient  features  of  this  technique. 

Because  of  its  geometrical  interpretation  the  problem  of  finding 
the  numerical  value  of  an  integral  of  a function  of  one  variable  is  often 
for  simplicity  called  quadrature.  One  method  of  quadrature  used  is  to 
evaluate  integrals  approximately  by  m:ans  of  a finite  number  of  values  of 
the  integrand.  In  many  cases  this  method  requires  less  work  tnan  other 
quadrature  methods. 

Quadrature  formulas  are  often  constructed  from  interpolating 
polynomials.  In  this  vay  we  can,  in  many  cases,  obtain  quadrature 
formulas  which  are  convenient  to  use  and  which  will  give  sufficiently 
accurate  results.  For  n arbitrary  points  Xj,x0,...,x  in  the  segment 
[a,b[,  we  can  construct  the  interpolating  polynomial  for  f(x)  by  the 
following  formulas 


f(x)  = P(x)  + r(x)  (A . 1 ) 

p(l°  ‘ j,  (»-^y(ic-(v  f(V  <A-2) 

W(x)  = (x-x,)  (x-x„)  ...  (x-x  ) (A. 3) 

i / n 

r(x)  = remainder  of  interpolation  (A. 4) 


[33]  V.  I.  Krylov,  "Approximate  Calculations  of  Integrals," 
Macmillan  Company,  New  York,  1962. 


-85- 


b 


The  exact  value  of  Che  Integral 


p(x)  f(x)dx  is 


b 

r 


p(x)f (x)dx 


a 

b b 

p(x)P(x)dx  + I 


p(x)r(x)dx 


(A.5) 


a a a 

where  p(x)  is  the  weight  function.  It  is  assumed  that  p(x)  is  a certain 
fixed  function,  which  is  measurable  on  [a,b]  and  is  not  Identically  the 
zero  function,  and  that  the  product  p(x)f(x)  is  integrable  on  [a,b].  If 
the  Interpolation  (A. 1)  is  sufficiently  accurate  so  that  the  remainder 
r(x)  is  small  throughout  the  interval  [a,b]  then  the  second  term  in  (A.5) 
can  be  neglected  and  the  approximate  equation  obtained  is 


b 

• n 

p(x)f  (x)dx  l A,  f (x.  ) 

k-1  K K 
a 


(A. 6) 


where 


\ 


b 

* 

pCx ) 

a 


W(x) 

(x-Xj^)  W'(xk) 


dx 


(A.  7) 


Quadrature  formulas  for  which  the  coefficients  have  the  form  (A. 7)  are 
called  interpolatory  quadrature  formulas.  Interpolatory  quadrature 
formulas  can  be  characterized  by  the  following  theorem. 


Theorem  1:  In  order  that  the  quadrature  formula  be  interpolatory  it  is 

necessary  and  sufficient  that  it  be  exact  for  all  possible  polynomials  of 
degree  n- 1 . 

b 


The  quadrature  formula 


n 

p(x)f(x)dx  ^ \ A^f(x^)  for  a fixed  n, 
k»l 


contains  2n  parameters  and  x^  (k-1 ,2, . . . ,n) . The  problem  is  to  select 
these  parameters  so  that  (A. 6)  will  be  exact  for  all  polynomials  of  the 
highest  degree  (l.e.  for  all  polynomials  if  degree  <_  k,  where  /.  is  as 
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large  as  possible).  The  choice  of  the  coefficients  for  any  arrange- 
ment of  x^  can  lead  to  an  equation  (A. 6)  which  is  exact  for  all  poly- 
nomials of  degree  <_  n-1.  This  requirement  completely  defines  the  coef- 
ficients A^:  (A. 6)  must  be  interpolatory  and  its  coefficients  must  be 
given  by  (A. 7) 

In  order  to  increase  the  precision  of  (A. 6)  the  choice  of  the 
points  is  still  at  our  disposal.  We  might  hope  that  for  some  choice 
of  these  points  the  degree  of  precision  can  be  Increased  by  n and  that 
the  formula  can  be  made  exact  for  all  polynomials  of  degree  2n-l. 

This  can  be  achieved  as  3hown  below.  The  conditions  which  must  be  satis- 

fied by  A^  and  x^  in  order  that  (A. 6)  will  be  exact  will  now  be  estab- 
1 ished. 

We  prefer  to  consider  the  polynomial  W(x)  « (x-x^) (x-x^) . . . (x-x^) 
instead  of  the  nodes  x^  themselves.  If  we  know  the  x^,  then  we  can  easily 
construct  the  polynomial  W(x).  Conversely,  if  we  know  the  polynomial 
W(x)  * x°  + a^x"  * + •••»  t^en  determining  the  roots  of  W(x)  will  give  us 

the  x^.  If  we  determine  W(x)  instead  of  x^  directly  then  we  must  be  care- 

ful that  the  roots  of  W(x)  will  be  real,  distinct  and  located  in  the  seg- 
ment [a,b]. 


Theorem  2 : If  (A. 6)  is  to  be  exact  for  all  polynomials  of  degree  (2n-l) , 

then  It  is  necessary  and  sufficient  that  (A. 6)  must  be  interpolatory  and 
that  the  polynomial  W(x)  be  orthogonal  with  respect  to  p(x)  to  all  poly- 
nomials Q(x)  of  degree  < n: 

b 

p(x)  W(x)  Q(x)dx  ■ 0 (A. 8) 

« 

a 
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4 


The  possibility  of  constructing  formulas  which  are  exact  for  all 
polynomials  of  degree  2n-l  is  related  to  the  existence  of  polynomials  W(x) 
of  degree  n which  possess  the  above  orthogonality  property.  If  the  weight 
function  p(x)  changes  sign  on  (a,b)  then  such  a polynomial  W(x)  may  not 
exist.  If  such  a polynomial  does  exist  its  roots  might  not  satisfy  the 
above  requirements.  Hence  it  will  be  assumed  that  the  weight  function 
p(x)  is  nonnegative  on  [a,b]. 

Theorem  3:  It  the  polynomial  P (x)  is  orthogonal  on  the  segment  [a,b] 

n 

to  all  polynomials  of  degree  less  than  n,  with  respect  to  gatWe 

weight  function  p(x),  then  all  the  roots  of  P^(x)  are  rea  «'inct 

and  lies  inside  [a,b]. 

Theorem  4:  If  p(x)  > 0 for  x > fa,b],  then  a quadrature  formula  (A. 6) 

which  is  exact  for  all  polynomials  of  degree  2n-l,  exists  for  all  n, 

and  cannot  he  exact  for  all  polynomials  of  degree  "'n. 

The  construction  of  quadrature  formulas  which  have  the  highest 

accuracy  is  now  discussed.  Consider  the  system  of  polynomials  P (x) ; 

n 

(n*l,2,...)  which  are  orthogonal  on  [a,b]  with  respect  to  the  weight 
function  p(x).  In  order  to  be  cefinite,  it  is  assumed  that  the  system 
is  normalized, 

b ( 0 for  m 4 n 

p (x)  P (x)  P (x)dx  - < (A. 9) 

n m ) 

a \ 1 for  m ■ n 

and  the  nth  degree  polynomial  of  an  orthonormal  system  can  be  assumed  to 
be  of  the  form 

P (x)  * a xn  + b x°  ^ + . . . 
n n n 
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The  nth  degree  polynomial  of  this  system  then  differs  from  W(x)  by  only 
a constant  multiple.  The  roots  of  P^Cx)  will  thus  be  the  nodes 
(k  ■ l,2,...n)  which  are  to  be  used  in  the  quadrature  formula.  The  coef- 
ficients A,  are  determined  by  (A. 7)  or  equivalently  by, 


\ 


b 


p(x) 


pn(x) 

O^Tp^OO  dx 


a 

In  order  to  calculate  A^  by  (A. 10)  the  Christoffel  Darboux  identity  is 
used.  Three  consecutive  polynomials  of  an  orthonormal  set  satisfy  a 
recursion  relation 


(A.  10) 


xP  (x)  * Pn+1  (x)  + P (x)  + P . (x) 

n a n+1  a a ,,  n a n-1 

n+1  n n+1  n 


(A. 11) 


from  which  the  Christof fel-Darboux  identity  can  be  deduced 


<’"t)  Jn  Pk<X)  Pk<t)  ‘ [pn+l<x)  Pn(t)  ‘ Pn(x)  Pn+l(t»  (A'12) 

k-0  n+1 


From  (A. 12)  it  can  be  written  that 

_ jn+1  1 _jn 1 

^ " *„  W P„+1<V  ' “n-1  W Pn-l(*t) 

Thus  from  the  second  equality  of  (A. 13)  it  is  clear  that  a quadrature 
formula  of  the  highest  accuracy  has  all  positive  coefficients. 

Next  the  magnitude  of  the  remainder  is  examined. 


(A. 13) 


Theorem  3:  If  f(x)  has  a continuous  derivative  of  order  2n  on  [a,b]  then 
there  exists  a point  $ in  [a,b]  for  which  the  remainder  of  the  quadrature 
formula  of  the  highest  accuracy  is 


rff) 


f(2n)(<fr) 
(2n) ! 


p(x)  (x)dx 


a 


(A. 14) 
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These  quadrature  formulas  of  highest  algebraic  accuracy  are  used  to 
evaluate  two  integrals  which  often  occur  in  the  specified  problem  of 
interest.  It  was  because  of  those  formulas  that  it  was  possible  to 
write  user-oriented  computer  programs  using  the  Sommerfeld  formulation 
that  do  not  require  unreasonable  amounts  of  time. 


a)  Constant  Weight  Function: 


The  formulas  of  Gauss  are  historically  the  first  formulas  of 
the  hi;  :.;st  algebraic  accuracy.  These  formulas  are  used  to  approximate 


the  integral 


f (x)dx 


where  [a,b]  is  a finite  segment  and  p(x)  F 1. 

By  a linear  transformation  we  can  transform  an  arbitrary  segment 
[a,b]  into  any  standard  segment  we  choose.  In  order  to  make  use  of  the 
symmetry  of  the  nodes  x^  and  coefficients  the  segment  would  be  made 
to  be  [-1,  +1].  Thus  it  will  be  assumed  that  (A. 15)  can  be  transformed 
into  the  form 


f -j^xjdx 


The  system  of  polynomials  which  are  orthogonal  on  [-1,  +1]  with  respect 
to  the  constant  weight  function  are  the  Legendre  polynomials, 


. ,n,  2 Nn 

P (X)  - -i-  d fe  -■1). 
n 2"nl  dx" 


The  quadrature  formula  of  the  highest  accuracy 


+1 

f(x)dx  = l A^n)  f(x^n)) 
-»  k.—  1 
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has  for  its  n nodes  the  roots  of  the  Legendre  polynomial  of  degree  n: 


’ 0 


(A.  19) 


The  coefficients  are  then  obtained  as 


Ai 


(n) 


[i  - (x<n))2]  [p;(*£°)]2 


(A. 20) 


jn  / | \ 2 

and  since  W(x)  * ; , P (x)  to  make  the  leading  term  of  W(x)  unity,  the 

f2n) ! n 

remainder  of  the  Gauss  formula  is  obtained  from  (A. 14)  as 

„2n+l  , . v 2 


r(f) 


(2n+l)  (2n) ! 1 (2n) ! 


■Mr)2 


(A. 21) 


where  $ is  a point  in  the  segment  [-1,  +1].  Values  of  x^  and  for 
different  n are  given  in  the  IBM-SSP  [34]. 


b)  Integrals  of  the  Form  J x exp(-x)  f(x)ds 

0 

The  system  of  polynomials  which  are  orthogonal  on  the  semi-infinite 
axis  0 £ x < “with  respect  to  the  function  xU  exp(-x)  is  the  system  of 
Chebyshev-Laguerre  polynomials 


,n 


L^(x)  * (-l)n  x “ exp(x)  (x0^  exp(-x)) 
n dx 


(A. 22) 


A quadrature  formula  of  the  highest  accuracy 


00 

t n 

xa  exp(-x)  f(x)dx  = J A^f(x^)  + r(f) 
a k=  1 


(A. 23) 


[34]  IBM  Scientific  Subroutine  Package  (360A-CM-03X)  Version  II, 
Programmer’s  Manual,  pp.  299-303. 


-91- 


must  have  as  its  nodes  the  roots  of  the  Laguerre  polynomial  (x) , and 

n 

the  coefficient  can  be  found  to  be 


(a  + n)  ! 

Values  of  and  for  a * 0 and  for  different  n are  given  in  the 
IBM-SSP  [35],  The  remainder  of  the  quadrature  formula  a * 0 is 


wiidre  $ is  a point  between 


r(f) 

[0, 


(n\)2 

(2n) ! 


f(2n) 


(<fr) 


°°] . 


(A.  24) 


(A.  25) 


[35]  IBM-SSP,  pp.  303-307. 
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Appendix  B:  ASYMPTOTIC  EVALUATION  OF  INTEGRALS  BY  THE  METHOD  OF 

STEEPEST  DESCENT 


The  method  of  steepest  descent  (or  the  saddle  point  method)  deals 
with  the  approximate  evaluation  of  integrals  of  the  form 


XSD^ 


F(£)  exp[-pf  (O  ]d£ 
C 


<B.1> 


for  large  values  of  p,  where  the  contour  C is  such  that  the  integrand 
goes  to  zero  at  the  ends  of  the  contour.  The  functions  f(£)  and  F(£) 
are  arbitrary  analytic  functions  of  the  complex  variable  £. 

The  basic  philosophy  of  the  method  of  steepest  descent  is  as 
follows:  A path  is  selected  in  the  complex  £ plane  in  such  a way  that 

the  entire  value  of  the  integral  is  determined  from  a comparatively 
short  portion  of  the  path.  Within  certain  limits,  the  contour  of  inte- 
gration C may  be  altered  to  such  a path  without  affecting  the  value  of  the 
integral.  Then,  the  integrand  is  replaced  by  another,  simpler  function, 
which  closely  approximates  the  integrand  over  the  essential  portions  of 
the  path.  The  behavior  of  the  new  integrand  outside  the  important 
portion  of  the  path  is  of  no  concern.  For  real  and  positive  value*  of  p 
and  for  a general  contour  C the  quantity  pf(C)  is  positive  on  some  parts 
of  the  path  and  there  are  other  regions  where  it  is  negative.  The  latter 
regions  are  more  important  since  the  integrand  is  larger,  and  in  these 
regions,  where  the  negative  of  the  Re[pf(E;)]  is  largest,  it  is  important 
to  reduce  oscillations.  A contour  is  chosen  along  which  the  imaginary 
part  of  [pf(C)l  is  constant  in  the  region  where  the  negative  of  its  real 
part  is  largest.  The  path  in  the  region  where  Re[pf(£)3  is  greatest  may 
be  chosen  so  that  Im[pf(01  varies  if  this  turns  out  to  be  necessary  to 
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complete  the  contour.  In  this  way,  the  oscillations  of  the  integral 
cause  the  least  trouble. 

This  general  philosophy  can  then  be  applied  to  evaluate  the 
integral  of  (B.l).  The  type  of  integral  encountered  in  this  problem 
is  of  the  form 


SD 


F(8)  exp[-  jkR  cos  ($-0)]dg 


(B.  2) 


where  F^  is  the  path  of  integration  as  shown  in  Fig.  5.  The  saddle 
point  occurs  at  f'(0  =0,  or  at  8 = 0 in  this  case.  The  path  of 
integration  is  now  determined  from 


cos(B-e)  = 1 - js  /k 


where  s is  real  and  -«•  <_  s + ®.  The  saddle  point  corresponds  to 
s = 0.  (B.  2)  can  now  be  expressed  in  the  following  form 


SD 


exp(-jkR) 


$(s)  exp(-R3  )ds 


(B.  3) 


where  now  9(b)  = F(8)  . The  details  of  evaluating  (B.3)  can  be  found 

in  [36]  and  only  the  result  is  quoted. 


rSD  = (^)1/2  exPU<*/4  - kR>iF(0) 

11  -■  {i  + + •••>  <B-4) 

Equation  (B.4)  is  inapplicable  if  there  is  a pole  near  the  saddle  point  0. 
However,  the  method  of  steepest  descent  can  be  modified  in  such  a way  that 

[36]  L.  B.  Felsen  and  N.  Marcuvitz,  "Radiation  and  Scatterin'  of  Waves," 
Prentice  Hall,  New  Jersey,  1973. 


I 
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the  presence  of  poles  is  taken  into  account  from  the  very  beginning, 
of  special  interest  in  the  analysis  will  be  an  integral  of  the  form 


I(kR)  = j F1(6)  exp[-jkR  cos(B-6)]d0  (B.5) 

F1 

where  is  a path  of  integration  in  the  complex  6 plane  as  shown  in 
Fig.  5 and  F^(8)  now  has  a pole  at  3p  which  is  near  0.  For  large  kR, 
the  following  method  is  applied:  The  pole  is  factored  out  from  F^(B) 

by  writing  F (8)  = . It  is  then  argued  by  Clemmow  [3/.  that 

i P-Pp 

sin  (— j- ) 

since  F(B)  has  no  singularities  in  the  vicinity  of  the  saddle  point  it 
may  be  removed  from  under  the  integral  sign  with  B equated  to  0.  Thu?, 
the  integral  I(kR)  can  be  written  as 


i.n(kR)  - F(to  f j.os(ere)l  dB 


T1  sin  (-^ — ) 


(B.  6) 


= F (0 ) 


exp[-ikr  cos  a]  , 

. ' ^-6P  d“ 
rQ  sin  ( — ^ ) 


(B.7) 


By  reversing  the  sign  of  a as 


ISD(kR)  = F(0) 


exp[-jkR  cos  a] 


0-a-Br 


da 


ro  8in  (— 


-) 


(B.  8) 


and  by  adding  (B.8)  and  (B.7)  and  then  dividing  by  two,  puts  (B.6) 
in  the  form 


[37]  P.  C.  Clemmow,  "The  Plane  Wave  Spectrum  Representation  of  Electro- 
magnetic Fields,"  Pergamon  Press,  New  York  1966,  pp.  46-58. 
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I (kR)  .2,1,1  F(0)  SSS.  °)  coS(f)da 

SD  2 J cos  a - cos  y 2 

ro 


where  y = 0 - Bp.  Changing  the  variable  of  integration  from  a to  t such 
that  t = /2  exp(-jn/4)  sin(y)  the  path  is  transformed  to  an  integral 
from  to  -H®.  Hence 

2 

r 

jb 


ISD(kR)  = 2h  exp 


[ j ( 3tt / 4 - kR)  ] F (0)  f g-x-Pir.k.RT  J 

j x + 1b 


where  b = sin  Application  of 


(B.9) 


exP(-XT  ) dT  = 

2 u.  IK2 
T + jb 


exp[  j (xb2 


- tt / 4)  ] erf  c (/j'xb2  ) 


from  Gradshteyn  and  Ryzhik  (p.  388,  3.466  No.  1)  to  (B.9),  yields 


I (kR)  = 2tt  exp[  j (n/2  - kR)  - W ] erfc  (jW) 


(B. 10) 


22  2 

where  W = - jkRb  = -j2kR  sin  ( — - — ).  Finally,  application  of  (B.4) 
to  the  infinite  integrals  encountered  in  the  Sommerfeld  formulation  yields 
the  reflection  coefficient  method  and  that  of  (B.10)  gives  very  good 
approximations  for  near  fields. 
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